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Abstract 


Currently,  the  Department  of  Defense  (DOD)  and  the  commercial  airline  industry 
are  utilizing  the  Instrument  Landing  System  (ILS)  during  aircraft  landings  for  precision 
approaches.  The  replacement  system  for  the  aging  ILS  was  thought  to  be  the  Microwave 
Landing  System  (MLS).  Instead,  use  of  the  Global  Positioning  System  (GPS)  is  now 
thought  to  be  a  viable  replacement  for  ILS  precision  approaches.  The  majority  of  current 
precision  landing  research  has  exploited  “stand-alone”  GPS  receiver  techniques.  This 
thesis  instead  explores  the  possibilities  of  using  an  extended  Kalman  filter  (EKF)  that 
integrates  an  Inertial  Navigation  System  (INS),  GPS,  Barometric  Altimeter,  Pseudolite 
and  Radar  Altimeter  for  aircraft  precision  approaches.  Thesis  results  show  that 
integrating  the  INS,  GPS,  Barometric  Altimeter  and  Radar  Altimeter  meets  Federal 
Aviation  Administration  (FAA)  requirements  for  a  Category  I  precision  approach.  Thesis 
results  also  show  that  integrating  the  INS,  GPS,  Barometric  Altimeter,  Radar  Altimeter 
and  a  single  Pseudolite  meets  FAA  requirements  for  a  Category  II  precision  approach. 
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AN  INTEGRATED  GPS/INS/BARO  AND  RADAR  ALTIMETER  SYSTEM 
FOR  AIRCRAFT  PRECISION  APPROACH  LANDINGS 


I.  Introduction 

Current  Department  of  Defense  (DOD)  and  commercial  aircraft  utilize  the 
Instrument  Landing  System  (ILS)  during  landings  for  precision  approaches  [19,22].  The 
replacement  system  for  this  aging  ILS  was  thought  to  be  the  Microwave  Landing  System 
(MLS).  Instead,  use  of  the  Global  Positioning  System  (GPS)  is  now  thought  to  be  a 
viable  replacement  for  aircraft  ILS  precision  approaches. 

The  DOD  and  Department  of  Transportation  (DOT)  are  under  extreme  pressure 
from  military  users  and  the  commercial  airline  market  to  "find  a  way"  to  certify  GPS  for 
precision  landing  approaches  [2,20,55,64,67].  Using  GPS  for  precision  landing 
approaches  instead  of  using  an  additional  MLS  receiver,  clearly  makes  sense  from  a 
cost/benefit  perspective.  First,  GPS  will  eventually  be  installed  on  all  the  DOD's  major 
weapon  systems  (aircraft,  ships,  etc.)  and  eventually  be  installed  on  commercial 
transports  [21].  Second,  if  a  MLS  receiver  does  not  have  to  be  purchased  by  the  DOD  or 
commercial  airlines,  billions  of  dollars  would  be  saved.  For  example,  if  only  the  F-16 
system  program  office  (SPO)  could  integrate  their  2000  aircraft  with  existing  GPS  and 
INS  properly,  the  F-16  SPO  could  conservatively  save  $250  million  dollars,  since  no 
additional  MLS  receiver  will  have  to  be  purchased  [37].  Studies  also  have  shown  the  use 
of  GPS  on  commercial  transports  will  save  airlines  an  estimated  6  billion  dollars 
annually.  The  DOD  would  also  benefit  from  substantial  dollar  savings.  Civil  aviation 
aircraft  can  also  benefit  from  the  use  of  GPS  and  GPS -aided  inertial  systems  for  non¬ 
precision  and  precision  approaches,  provided  that  the  total  system  cost  is  not  prohibitive 
for  such  aircraft.  In  an  effort  to  reduce  costs,  airlines  are  looking  for  ways  GPS  can  save 
operating  costs.  The  author  believes  that  in  order  for  GPS  to  be  certified  on  civil  and 
DOD  aircraft,  not  only  must  GPS  meet  the  accuracy  requirements,  but  the  integrity  of  the 
GPS  must  not  be  in  question.  The  author  also  believes  that  a  GPS  combined  with  an  INS 
can  meet  both  the  accuracy  and  integrity  requirements  of  the  FAA  for  precision 
approaches. 

If  one  currently  looks  at  how  the  majority  of  FAA  certified  precision  approach 
aircraft  are  equipped  with  avionics  black  boxes,  one  would  typically  see  the  following 
shopping  list  of  Table  1-1: 
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Quantity  (each) 

Nomenclature 

1 

Barometric  altimeter 

1  to  3 

Inertial  Navigation  System 

1 

Radar  altimeter 

1* 

Global  Positioning  System* 

— 

*  =  predicted  to  be  installed  by  1998 

Table  1  -  1 .  Typical  avionics  equipment  found  on  commercial  or  DOD  aircraft 

Present-day  INS  systems  on  board  most  commercial  and  DOD  aircraft  are  medium 
accuracy  (1.0  nm/hr  CEP)  stand-alone  systems  at  a  cost  of  about  $100,000  per  unit.  The 
author  believes  that  the  present-day  stand-alone  INS  could  be  replaced  by  a  single 
integrated  INS-GPS  black-box  which  will  utilize  barometric  and  radar  altimeter 
measurements  to  meet  and  exceed  FAA  requirements  at  a  reduced  cost;  see  Table  1-2 
for  cost  comparisons  [37,62,61].  The  reduction  in  cost  is  found  by  replacing  the  current 
medium  1.0  nm/hr  accuracy  INS  with  a  less  accurate,  2-4  nm/hr  INS  tightly  coupled  with 
a  carrier-phase  differential  GPS.  Table  1-2  does  not  show  additional  integration  cost 
savings  of  an  off-the-shelf  GPS -INS  which  will  already  contain  a  host  Kalman  filter 
algorithm.  A  further  advantage  to  such  a  proposed  system  is  that  the  INS  can  reduce 
carrier-phase  ambiguities,  and  GPS  integrity  monitoring  techniques  can  be  employed 


Black-box 

Single  INS 

Triple  INS 

Single 

Embedded  INS- 

GPS 

Triple 

Embedded  INS- 

GPS 

Baro  altimeter 

$20K 

$20K 

$20K 

$20K 

INS 

$100K 

$300K 

— 

— 

Radar  altimeter 

$30K 

$30K 

$30K 

$30K 

GPS 

$25K 

$25K 

— 

— 

Embedded  INS- 

GPS 

— 

— 

$60K 

$180K 

Total 

$175K 

$375K 

$110K 

$230K 

Table  1-2.  Cost  Comparison  of  Stand-Alone  vs.  Embedded  INS-GPS 
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using  a  centralized  Kalman  filter.  Another  advantage  of  an  embedded  INS-GPS  black¬ 
box  is  that  it  contains  a  Kalman  filter  "tuned"  to  the  specific  INS-GPS  hardware. 

In  order  for  GPS  to  replace  ILS  for  non-precision  and  precision  approaches,  areas 
associated  with  accuracy,  coverage,  integrity,  availability  and  aircraft  integration  must  be 
studied  and  understood.  Because  a  GPS  differential-carrier-phase  model  was  unavailable 
at  the  onset  of  this  thesis,  the  author  will  instead  use  a  P-Code  GPS  model,  with  and 
without  the  accurate  radar  altimeter.  The  beauty  of  the  radar  altimeter  is  it  is  designed  to 
do  one  thing  -  -  give  height  measurements  referenced  above  the  ground.  It  is  no  secret 
that  during  an  aircraft  landing,  incorporation  of  height  above  ground  data  is  of  utmost 
importance,  yet  most  commercial  and  military  aircraft  do  not  integrate  radar  altimeter 
information  using  a  Kalman  filter.  Instead,  the  author  is  aware  that  most  military  and 
commercial  airlines  have  purchased  the  radar  altimeter,  but  it  is  only  used  in  a  stand¬ 
alone  mode.  Another  reason  to  integrate  the  radar  altimeter  height  above  ground 
measurements,  is  that  typically  code  phase  GPS  suffers  most  in  the  vertical  channel  due 
to  satellite  vehicle  geometry  than  when  compared  with  GPS  horizontal  position  errors. 
The  author  is  aware  that  not  all  aircraft  (particularly  civil  aviation  aircraft)  can  afford  the 
use  of  a  radar  altimeter;  instead,  the  use  of  a  GPS  pseudolite  (used  during  aircraft  final 
approach)  will  be  compared.  This  thesis  will  take  a  first-cut  look  at  comparing  the 
following  integrated  systems  containing  a  low  and  medium  accuracy  INS: 


Case  I 

Case  11 

Case  III 

Case  IV 

CaseV 

Case  VI 

CaseVII 

CaseVIII 

Case  IX 

Barometric 

Barometric 

Barometric 

Barometric 

Barometric 

Barometric 

Barometric 

Barometric 

Barometric 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

Altimeter 

0.4  nm/hr 

0.4  nm/hr 

2.0  nm/hr 

2.0  nm/hr 

4.0  nm/hr 

4.0  nm/hr 

0.4  nm/hr 

2.0  nm/hr 

4.0  nm/hr 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

CEP  INS 

P-Code 

P-Code 

P-Code 

P-Code 

P-Code 

P-Code 

P-Code 

P-Code 

P-Code 

GPS 

GPS 

GPS 

GPS 

GPS 

GPS 

GPS 

GPS 

GPS 

Radar 

Radar 

— 

Radar 

Pseudolite 

Pseudolite 

Pseudolite 

Altimeter 

Altimeter 

Altimeter 

Table  1-3.  Case  I-IX  Thesis  Integration  Comparisons 
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1.1 


Background 


After  much  review  of  the  past  years'  technical  publications,  it  appears  that  very 
few  "integrated"  GPS/INS  studies  have  been  accomplished  with  regard  to  implementing 
precision  landing  approaches.  Instead,  the  majority  of  research  has  exploited  "stand¬ 
alone"  GPS  receiver  techniques  (Differential  GPS  (DGPS)  or  Carrier-phase  GPS 
(CPGPS))  [53,66,26,5,51,12,30,34].  Only  two  publications  described  an  integrated 
mechanization  using  an  Inertial  Navigation  System  (INS)  [30,34]. 

No  technical  papers  were  noted  using  a  radar  altimeter,  INS  and  GPS 
measurements  with  a  Kalman  filter.  The  author  chose  to  use  a  radar  altimeter  because  the 
radar  altimeter  can  provide  vertical  (altitude)  measurements  to  the  Kalman  filter  which 
otherwise  would  not  be  utilized.  Integration  techniques  for  precision  approaches 
normally  fail  because  the  vertical  channel  accuracy  requirements  are  not  met  (despite  the 
use  of  a  barometric  altimeter  to  stabilize  the  vertical  channel  of  the  INS).  The  accuracy  of 
most  radar  altimeters  can  meet  the  vertical  precision  approach  requirements  typically  to 
Category  I  requirements  (see  Section  1.2,  Table  1  -  4)  [28].  Another  benefit  of  the  radar 
altimeter  is  it  can  be  used  in  ground  collision  avoidance  algorithms.  Because  the  radar 
altimeter  is  equipped  on  most  commercial  airline  and  DOD  aircraft,  and  because  the  radar 
altimeter  shows  great  potential  for  aiding  an  aircraft  during  precision  approaches  using  a 
Kalman  filter,  it  was  chosen  to  be  utilized  in  this  thesis.  The  author  believes  that 
pseudolites  (active  pseudo-"satellites"  precisely  stationed  along  the  approach  runway 
path)  show  great  potential  in  allowing  the  host  aircraft's  Kalman  filter  access  to  a  superb 
ranging  device  to  eliminate  vertical  errors.  These  are  currently  being  investigated  by  a 
number  of  researchers  [12],  and  so  they  will  also  be  considered  in  this  performance 
analysis.  It  was  also  noted  that  though  a  wide  variety  of  integration  techniques  exist 
(fuzzy  logic,  neural  networks,  least  squares,  etc.),  only  Kalman  filter  integration 
techniques  were  found  to  be  actively  in  use  for  current  aircraft  precision  approach 
research  and  development. 

1.2  Key  Terms 

Global  Positioning  System:  A  satellite-based  navigation  and  time  system  designed, 
developed  and  maintained  by  the  U.S.  Department  of  Defense  (DOD)  [17].  Accuracy: 
Stand-alone  (no-aiding)  GPS  is  typically  16  meters  (military)  and  100  meters 
(commercial)  spherical  error  probable  (SEP). 
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Differential  GPS  (DGPS):  A  ground-based  GPS  receiver  (accurately  surveyed)  uplinks 
error  corrections  to  nearby  (<  150  nautical  miles)  aircraft,  thus  reducing  common  errors. 
Accuracy:  3  meters  (military  and  commercial)  SEP. 

Carrier- Phase  GPS  (CPGPS):  A  new  receiver  technique  which  is  able  to  measure  the 
incoming  satellite-transmitted  GPS  signal  to  a  fraction  of  a  wavelength.  Accuracy:  <  30 
centimeters  (military  and  commercial),  1-a. 

Inertial  Navigation  System  (INS):  A  self-contained  dead-reckoning  system  that  utilizes 
internal  gyroscopes  and  accelerometers  to  navigate.  Typical  medium  accuracy:  1-3 
nautical-miles/hour  circular  error  probable  (CEP). 

Kalman  Filter:  A  recursive  computer  algorithm  that  uses  sampled-data  measurements  to 
produce  optimal  estimates  of  states  of  a  dynamic  system,  under  the  assumptions  of  linear 
system  models  and  white  Gaussian  noise  models.  Developed  by  R.E.  Kalman  in  the  early 
1960s,  A  New  Approach  to  Linear  Filtering  and  Prediction  Problems  [35], 

Aircraft  Precision  Landing:  Formally  defined  by  the  Federal  Aviation  Administration 
(FAA)  in  [22]  as  Category  I,  II  or  HI  precision  approach.  See  Table  1  -  4  [36]. 


Precision  Approach  Parameters 

Category 

1 

(in  feet,  all  1 -sigma  values) 

Azimuth 

+/-  28.1 

Elevation 

+/-  6.8 

II 

+/-  8.6 

+/-  2.8 

III 

+/-  6.8 

+/- 1.0 

_ 1 

Table  1-4.  Precision  Approach  Accuracy  Requirements  at  Decision  Heights 

Instrument  Landing  System  (ILS):  Current  land-based  navigation  aide  used  to  guide 
aircraft  safely  for  final  approach  airport  landings  [22,36]. 

Microwave  Landing  System  (MLS):  Proposed  land-based  replacement  navigation  aide 
for  the  ILS  [36,43], 
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1.3  Literature  Review 


The  following  is  a  brief  discussion  of  eight  technical  papers  published  in  1993- 
1994  [53,66,26,5,51,12,30,34].  The  eight  technical  papers  were  reviewed  and 
summarized  in  three  categories:  (1)  DGPS  (stand-alone),  (2)  CPGPS  (stand-alone)  and 
(3)  GPS,  DGPS  or  CPGPS  integrated  with  an  INS  using  a  Kalman  filter. 

1.3.1  DGPS  ( Stand-Alone ) 

[53].  Commercial  code  tracking  differential  GPS  (DGPS)  landing  system,  using 
narrow  correlator  receivers  in  the  ground  reference  station  and  avionics,  configured  to 
drive  an  ILS  autoland  flight  control  system  with  ILS  "look  alike"  deviation  signals, 
successfully  guided  a  BOEING  737  to  31  successful  "hands  off'  landings.  No  landings 
are  aborted  because  of  equipment  failure,  and  conservative  estimates  of  lateral  and 
vertical  total  system  error  fell  within  Category  III  for  both  the  approach  and  touchdown 
segments  of  the  landings  with  substantial  margin.  Pilots  commented  that  DGPS  approach 
paths  seemed  noticeably  straighter  than  what  they  had  experienced  with  ILS  coupled 
approaches.  Performance  results  were  based  on  predicted  FAA  "tunnel"  requirements. 

[66].  A  commercial  code  tracking  Novatel  Model  95 1R  GPS  card  is  utilized, 
along  with  DGPS  corrections.  A  total  of  43  approaches  are  flown  at  Atlantic  City,  NJ. 
The  Novatel  receiver  is  a  10-channel,  narrow  correlator  which  outputs  data  at  a  5  Hz  rate. 
Category  I  vertical  requirements  are  met. 

1.3.2  CPGPS  ( Stand-Alone ) 

[26].  The  theoretical  use  of  carrier-phase  for  Category  I,  II  and  in  approaches  is 
discussed.  Mathematical  equations  describing  single,  double  and  triple  differencing 
techniques  are  derived  for  use  with  pure  carrier-phase  GPS  receivers. 

[5].  The  theoretical  use  of  carrier-phase  GPS  and  solving  the  carrier-phase  cycle 
ambiguity  is  discussed.  The  author  describes  situations  in  which  use  of  only  LI  band 
cannot  guarantee  Category  HI  integrity  and  reliability.  Use  of  pseudolites  (GPS 
"satellites"  based  on  the  ground  at  precisely  known  locations)  to  carry  out  a  continuous 
landing  approach  is  mentioned. 

[51].  The  NASA  Ames  Research  Center  conducted  theoretical  research  to 
develop  and  demonstrate  carrier-phase  DGPS  algorithms  for  approach  and  landing.  The 
theoretical  research  is  put  to  test  by  actually  flight  testing  the  algorithms  and  navigation 
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systems  using  a  King  Air  200  aircraft  and  an  Ashtech  P-12  GPS  receiver.  In  1 1  of  12 
approaches,  integer  ambiguities  (See  Note  below)  are  resolved  at  a  minimum  distance  of 
2.7  km  from  landing,  with  the  solution  being  maintained  through  touchdown.  However, 
in  5  of  the  12  final  approach  and  landing  segments,  wrong  integers  are  determined. 
Continuing  work  is  necessary  to  achieve  "robust"  integer- ambiguity  resolution  without 
the  use  of  pseudolites. 

Note:  The  GPS  carrier  frequencies  (LI  and  L2)  are  in  the  L-band  (1.2  to  1.6  GHz)  [27]. 

The  wavelength  range  of  the  L-band  is  therefore  18  to  25  centimeters.  The  observed 
carrier-phase  quantity  is  a  measure  of  the  phase-shift;  that  is,  it  represents  only  a 
fractional  part  of  one  wavelength  [27].  The  total  phase-range  measurement  at  some  time 
epoch  t  would  then  be  equal  to  this  fractional  part,  <3>frac,  plus  an  integer  number  of  phase 
cycles  from  the  initial  time,  to,  to  the  time  epoch  t  (which  is  continuously  measured  and 
compensated  by  the  receiver),  Oint(to,t),  plus  an  integer  phase  ambiguity  term  N  [27], 

Note  that  <bmt(t0,t0)  represents  the  receiver's  initial  estimate  of  the  number  of  integer  phase 
cycles  [27].  The  integer  ambiguity  term,  sometimes  referred  to  as  the  cycle  ambiguity,  is 
the  difference  between  the  true  integer  count  at  time  t  and  the  current  integer  count  at 
time  t  measured  by  the  receiver  [27].  The  total  phase  can  then  be  represented  from  [27]: 

^totalW  =  ^fracW  +  ^intC^o.O  +  N(f)  (1.1) 

For  more  information  about  carrier-phase  GPS  terminology,  see  [27,6], 

[12],  Stanford  University  provided  flight  test  results  of  a  carrier-phase  GPS,  using 
pseudolites  along  the  airport  approach  path.  Flight  test  results  are  compared  to  data  from 
independent  laser  measuring  equipment.  The  "Kinematic  GPS"  achieved  centimeter 
accuracy  throughout  all  flight  tests  (between  5  and  30  centimeters,  1-a  =  one  standard 
deviation).  The  estimated  cost  to  outfit  a  runway  with  pseudolites  for  Category  ID 
precision  is  stated  at  less  than  $100K. 

1.3.3  GPS ,  DGPS  or  CPGPS  Integrated  with  an  INS  Using  a  Kalman  Filter 

[30].  A  carrier-phase,  double  differential  GPS  and  an  INS  are  integrated  using  a 
centralized  Kalman  filter  (post-processed  after  flight).  Through  the  analysis  of  flight  test 
data,  it  is  shown  that  the  GPS/INS  position  data  agrees  at  the  centimeter-level  (1-a)  with 
the  GPS-only  carrier-phase  receiver  position.  The  ability  of  the  INS  to  detect  and  correct 
cycle  slips  is  also  demonstrated  for  system  integrity  and  reliability  concerns. 
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[34].  An  INS  and  a  commercial  GPS  are  integrated  using  a  Kalman  filter  and 
flight  tested  in  Germany.  The  requirements  for  automatic  landing  systems  (accuracy, 
integrity  and  availability)  are  discussed  in  this  paper.  INS/GPS  integration  performance 
results  are  on  the  order  of  1  meter  accuracy.  The  performance  of  this  system  would  meet 
Category  III  ELS  specifications  and  could  be  used  at  airports  not  equipped  with  a  ground- 
based  GPS  or  ILS. 

1.4  Problem: 

In  the  case  of  aircraft  integration,  very  few  "integrated"  GPS/INS  studies  have 
been  accomplished  with  respect  to  implementing  precision  landing  approaches.  Instead, 
the  majority  of  research  has  explored  stand-alone  GPS  receiver  (mostly  using  uplinked 
GPS  differential  signals  to  reduce  atmospheric  and  selective  availability  errors).  The 
benefit  of  "stand-alone"  is  most  definitely  the  cost-savings  of  such  a  scheme.  The 
disadvantage  is,  GPS  as  a  stand-alone  is  not  an  error-proof  system.  It  is  a  complicated 
system  devised  of  three  main  segments:  (1)  space,  (2)  control  and  (3)  user.  Of  the  three 
parts  the  “user”  segment  is  most  vulnerable  and  is  susceptible  to  using  bad  satellite 
vehicle  (SV)  range  data  or  terrorist  jamming  (or  unintentional  jamming)  if  it  is  operating 
as  a  "stand-alone"  system.  During  a  precision  landing,  the  navigation  solution,  that  is,  the 
input  control  signal(s)  to  the  autopilot  (for  autoland),  must  be  of  highest  integrity  to  be 
FAA  certified. 

Instead  of  stand-alone  GPS  techniques,  this  thesis  will  integrate  stand-alone 
systems  (ENS/GPS/Baro/Radar  Altimeter)  using  an  extended  Kalman  filter.  With  this 
powerful  integration  of  available  sensor  measurements,  it  then  is  possible  to  perform 
residual  monitoring  [40]  (a  possible  follow-on  to  this  thesis).  Residual  monitoring  would 
be  utilized  as  a  “reasonableness”  check  of  the  incoming  GPS  measurements  so  that  the 
overall  navigation  solution  processed  by  the  INS/GPS  Kalman  filter  integration  is  of 
highest  integrity.  In  short,  use  of  an  INS/GPS  Kalman  Filter  integration  and  residual 
monitoring  would  be  exploited  to  try  to  meet  all  FAA  requirements  with  regard  to 
accuracy,  integrity,  coverage,  availability  and  continuity  for  precision  landing  operations. 
This  thesis  will  only  look  at  using  an  extended  Kalman  filter  to  meet  FAA  precision 
approach  accuracy  requirements. 

The  user's  GPS  receiver  would  make  use  of  the  good  high-frequency  INS 
information.  Also,  the  highly  accurate  low-frequency  GPS  information  would  be  used  to 
offset  the  long  term  errors  of  the  INS  (i.e.,  errors  due  to  gyro  drift).  A  GPS  and  an  INS, 
when  integrated  properly,  are  nearly  perfect  complements  of  one  another. 
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1.5  Problem  Statement: 


This  thesis  concentrates  on  setting  up  reliable  models  for  a  medium  (1  nm/hr 
CEP)  and  low  accuracy  INS  (4  nm/hr  CEP),  GPS  (4  channel  receiver,  5m  vertical  and  4m 
horizontal  precision,  la),  Baro  Altimeter  (50  -  150ft,  la),  and  Radar  altimeter  (1%  of 
altitude  ±lft,  la).  This  thesis  also  develops  a  generic  precision  approach  flight  profile 
(using  PROFGEN  [47])  that  encompasses  a  majority  of  aircraft  types.  Lastly,  this  thesis 
utilizes  a  single  ground-based  SV  (pseudolite)  and  available  true  post-processed 
ephemeris  data  [14],  instead  of  prior  FORTRAN  orbit  functions  used  at  AFIT 
[58,3,45,50,27].  Once  all  the  above  elements  are  in  place,  the  Multimode  Simulation  for 
Optimal  Filter  Evaluation  (MSOFE)  [48]  is  utilized  to  perform  extended  Kalman  Filter 
integration  analysis. 

1.6  Past  Research: 

The  past  research  at  the  Air  Force  Institute  of  Technology  (AFIT)  began  with  the 
generation  of  computer  models  for  INS,  GPS  and  Range/Range-Rate  System  (RRS),  as 
well  as  the  development  of  an  integration  scheme  to  blend  the  information  from  these 
three  independent  sources  into  a  single  navigation  solution.  The  choice  for  blending 
measurement  information  is  an  extended  Kalman  filter  (EKF).  The  overall  navigation 
system  developed  by  the  early  research  became  known  as  the  Navigation  Reference 
System  (NRS)  [45,50]  model.  This  thesis  will  use  pieces  of  the  NRS  model  previously 
developed  by  researchers  at  AFIT.  A  fighter  profile  provided  by  Wright  Laboratories 
[47]  is  used  initially  to  ensure  that  the  modified  NRS  model,  (from  this  point  on,  we  will 
call  the  modified  NRS  model  the  “Landing  System  Model”  (LSM))  closely  resembles 
prior  AFIT  research  results. 

1.7  Scope: 

Itemizing  the  majority  of  the  thesis  tasks  yields: 

a.  Review  prior  AFIT  thesis  of  Mosle  (GE  93D)  and  Negast  (GE  91D). 

b.  Research  current  aircraft  ELS  Category  I,  II  and  IH  precision  approach  techniques  and 
performance  specifications.  (Translate  required  performance  specifications/capabilities 
into  the  development  of  the  INS/GPS  integration  for  this  project). 
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c.  Interview  pilots  and  Air  Traffic  Control  personnel  who  have  "real-world"  experience 
in  precision  landing  modes  of  operation.  [Literature  readings  can  never  give  all  the 
information  needed  to  understand  a  problem  fully.  Since  real-world  knowledge  (from 
military  and  commercial  pilots)  is  available,  it  makes  sense  to  search  it  out]. 

d.  Create  a  realistic  scenario  and  flight  profiles  of  an  autoland  at  Wright-Patterson  AFB, 
OH.  PROFGEN  [47]  will  be  utilized  to  create  a  trajectory  that  will  represent  a  generic 
transport  aircraft  from  take-off  to  a  precision  approach  landing.  This  flight  profile  will  be 
called  “Tanker”  since  its  attributes  closely  resemble  those  of  a  KC-135  tanker  (Boeing 
707)  aircraft. 

e.  Find/synthesize  a  "truth  model”  (a  complete,  complex  mathematical  error  state  model 
that  portrays  true  system  behavior  very  accurately).  An  error  state  model  is  chosen,  as 
opposed  to  a  total  state  model,  because  an  error  model  is  more  adequately  represented  as 
linear  and  because  it  involves  slower  dynamics  than  a  total  state  model  for  an  INS; 
essentially  all  terrestrial  navigation  system  Kalman  filters  are  based  upon  error  state 
models  [40].  The  NRS  model  of  Mosle  is  used  by  extracting  from  the  NRS  truth  and 
filter  models,  the  doppler-aided  velocity  measurements  and  the  transponder  (RRS) 
measurements.  Replace  the  current  MSOFE  NRS  algorithms  that  "simulate"  satellite 
vehicle  (SV)  ephemeris  data  by  instead  using  true  SV  ephemeris  data  acquired  directly 
from  the  National  Geodetic  Office.  Modify  Mosle’ s  MSOFE  NRS  code  to  accept  radar 
altimeter  measurements.  This  modified  NRS  code  that  accepts  radar  altimeter 
measurement  will  be  called  the  "Landing  System  Model  (LSM).  (Unfortunately,  due  to 
the  many  unknowns  of  the  proposed  integration  mechanization  and  the  time  limitations  of 
this  thesis,  only  a  first-cut  approach  at  a  radar  altimeter  model  is  incorporated). 

f.  Linearize  the  truth  model  (if  necessary). 

g.  Conduct  a  performance  analysis  of  each  proposed  Kalman  filter  (of  item  (e)  above), 
driven  by  measurements  derived  from  the  truth  model  of  the  real  system.  [Repeat  step  (h) 
below  as  necessary]. 

h.  "Tune"  each  filter  (of  step  (g)  above)  to  provide  the  best  possible  performance. 

i.  Generate  a  Monte  Carlo  analysis  of  designs  (of  step  (h))  that  show  the  most  promise. 

j.  Conduct  a  performance/computer  loading  tradeoff  analysis  -  -  then  select  a  design. 
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k.  Implement  the  chosen  design  to  be  used  in  the  final  system  against  a  realistic  profile  of 
step  (d).  Perform  checkout  and  final  tuning  of  the  system. 


l.  Analyze  results  and  identify  real  world  implementation  issues. 

m.  Compare  performance  results  found  from  this  thesis  to  benchmark(s)  of  today's  ILS, 
MLS,  and  GPS  (stand-alone). 

1.8  Assumptions: 

All  theses  are  limited  by  the  assumptions  made,  and  no  research  can  be  adequately 
evaluated  unless  these  assumptions  are  clearly  defined  [45].  This  section  outlines  the 
assumptions  that  have  been  made  in  this  thesis. 

1.  All  work  has  been  conducted  through  computer  simulation.  The  "real"  world  in  the 
simulation  is  modeled  as  a  full-order  truth  error-state  model.  The  full-order  truth  and 
filter  models  are  presented  Chapter  3. 

2.  The  INS  platform  is  assumed  to  be  stabilized  with  a  barometric  (baro)  altimeter.  An 
INS  platform  is  unstable  without  an  outside  measurement  source  in  the  vertical  channel 
[9].  While  a  baro  altimeter  is  not  the  only  way  to  stabilize  a  platform,  it  is  a  commonly 
used  method.  The  use  of  the  baro  altimeter  is  included  in  the  modeling  of  the  system. 

The  majority  of  commercial  and  military  aircraft  utilize  a  radar  altimeter  in  a  stand-alone 
mode  for  terminal  approaches  to  a  runway.  This  thesis  will  instead  exploit  the  radar 
altimeter  as  an  independent  measurement  device  feeding  an  extended  Kalman  filter.  The 
radar  altimeter  measurements  will  be  utilized  at  altitudes  below  3000  feet  above  ground 
level  (AGL).  In  summary,  this  thesis  will  use  both  the  barometric  and  radar  altimeter 
measurements. 

3.  A  sample  period  of  one  second  has  been  chosen  (unless  otherwise  noted)  for  the  EKF. 
The  sample  period  refers  to  how  often  the  GPS  and  radar  altimeter  measurements  will  be 
brought  into  the  EKF.  Past  AFIT  research  has  used  a  variety  of  sample  periods,  varying 
from  two  to  ten  seconds  [45,50],  The  decision  to  use  one  second  sample  period  is  based 
primarily  on  the  typical  availability  of  the  GPS  measurement  in  the  real  world.  Though 
the  author  is  aware  of  a  few  GPS  receivers  which  output  measurements  at  a  rate  of  ten 
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times  a  second  (10  Hz),  a  one  second  sample  period  is  chosen  as  a  good,  representative 
design  choice. 

4.  The  computer  simulations  have  been  developed  using  a  program  called  Multi-mode 
Simulation  for  Optimal  Filter  Evaluation  (MSOFE)  [48].  MSOFE  is  well-established  Air 
Force  software  to  develop  and  test  Kalman  filter  algorithms. 

5.  The  computer-simulated  flight  profile  has  been  generated  by  the  program  PROFGEN 
[47],  PROFGEN  is  designed  to  work  with  MSOFE  to  provide  the  necessary  data  files  to 
simulate  dynamic  flight  profiles. 

6.  The  plotted  outputs  are  generated  by  the  commercial  software  package  MATLAB 
[60]. 

7.  The  SV  ephemeris  data  using  System  Effectiveness  Model  (SEM)  [18]  software  was 
obtained  from  the  Coast  Guard  BBS.  The  ephemeris  data  is  post-processed  by  the  U.S. 
Department  of  Commerce,  National  Geodetic  Information  Branch  [14]. 

8.  Ephemeris  data  was  incorporated  into  PROFGEN’ s  binary  output  “FLIGHT”  profile 
by  making  adaptations  to  existing  [57]  FORTRAN  source  code. 

9.  The  four  SVs  chosen  to  range  during  operation  of  MSOFE  and  the  FLIGHT  profile  are 
chosen  based  on  the  indicated  results  of  the  System  Effectiveness  Model  (SEM3.6) 
software  from  [  1 8]  based  on  position  dilution  of  position  (PDOP)  criteria  less  than  two 
(2). 

10.  The  simulation  software,  MSOFE  and  MATLAB,  has  been  coded  to  run  in  double 
precision  to  increase  the  numerical  stability  and  precision  of  the  simulation.  MSOFE 
software  utilizes  a  U-D  factorization  algorithm  to  increase  the  numerical  stability  in  the 
Kalman  filter  measurement  update  equations  [48,41]. 

11.  The  MSOFE  runs  are  conducted  using  15-run  Monte  Carlo  analyses.  While  a  larger 
batch  size  for  the  Monte  Carlo  analysis  would  be  preferable,  this  value  has  been  chosen  to 
keep  the  computational  burden  of  the  thesis  within  reasonable  bounds,  while  maintaining 
adequate  confidence  that  the  resulting  sample  statistics  properly  reflect  the  true 
underlying  statistics. 


1  - 12 


12.  Taylor  series  approximations  truncated  at  first  order  are  used  for  linearizing 
nonlinear  equations  in  the  NRS  and  LSM  filter.  Perturbations  about  a  nominal  trajectory 
will  be  established  in  each  case. 

13.  It  will  be  assumed  for  this  thesis  that,  when  radar  altimeter  measurements  are 
available,  the  earth’s  surface  will  be  modeled  as  flat  and  referenced  approximately  to  the 
INS  indicated  altitude  (referenced  to  WGS-84  ellipsoid).  This  assumption  will  definitely 
have  to  be  “upgraded”  to  a  more  realistic  radar  altimeter  scenario  at  a  later  time  by 
possibly  using  a  database  that  contains  “height  of  terrain”  for  specific  locations  on  the 
earth. 

14.  The  INS  will  have  had  a  “normal”  8-minute  alignment  and  nominal  flight  of  sixty 
(60)  minute  duration  prior  to  the  terminal  approach  phase  under  investigation.  “Normal” 
also  means  the  INS  has  not  been  degraded  nor  enhanced  by  any  means. 

15.  Four  SV  are  always  available,  with  an  average  PDOP  of  2.1. 

16.  The  transport  aircraft  flight  profile  will: 

a.  Always  be  at  less  than  0.9g  during  entire  flight. 

b.  Have  a  takeoff  speed  of  150  knots. 

c.  Have  a  landing  speed  of  133  knots  at  a  3  degree  glideslope. 

d.  Airspeed  above  10,000  will  always  be  greater  than  250  knots. 

e.  Change  altitude  at  a  rate  4000  ft/min  (maximum). 

f.  Follow  the  approach  plate  of  Figure  1  -  1  [15]. 

The  approach  plate  of  Figure  1  -  1  is  rather  cryptic  unless  one  is  a  trained  pilot.  Figure  1  - 
1  shows  two  top  views  of  the  Wright-Patterson  runway,  and  a  side  view  for  a  landing 
specifically  on  runway  23R.  The  aircraft  is  transitioned  to  the  Wright-Patterson  (FFO) 
procedure  track,  approximately  1 1  nm  from  the  runway.  The  aircraft  then  follows  the 
arc-turn  to  heading  233°.  During  inbound  transition,  the  aircraft  must  maintain  a 
minimum  altitude  of  3000  ft  MSL.  At  approximately  5.6  nm  from  the  VORTAC  DME 
station,  the  aircraft  maintains  a  minimum  of  2200  ft  MSL  which  is  the  glide  slope 
intercept  altitude.  At  approximately  0.5  nm  from  the  runway  end,  the  pilot's  VORTAC 
DME  would  indicate  2.1  nm.  The  aircraft  continues  its  3°  ILS  glideslope  to  touchdown. 
Figure  1  -  1  shows  that  FFO  is  only  certified  for  a  Category  I  precision  approach.  The 
lower  left  comer  of  Figure  1-1  shows,  "S-ILS  23R  1025/24  200."  This  information 
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states  that  it  is  a  "Precision,  straight-in  to  Runway  23  (right-hand  side);  the  decision 
height  (DH)  MSL  is  1025  ft  with  prevailing  visibility  (runway  visible  range  in  100's  feet) 
of  24."  "Height  of  DH  above  touchdown  zone  (HAT)  is  200  feet."  For  a  precision  ILS 
approach,  the  pilot: 

•  Transitions  to  the  ELS  Localizer  Course  from  the  published  approach  procedure 

•  Tunes  the  ILS  and  monitors  the  proper  identifier  during  the  entire  approach 

•  Sets  the  published  localizer  course  prior  to  localizer  course  interception 

•  Accomplishes  the  Approach 

•  Once  the  localizer  course  is  intercepted,  maintains  glide  slope  interception 
altitude  until  reaching  the  glide  slope  intercept  point. 

•  Maintains  a  complete  instrument  cross-check  throughout  the  approach,  with 
increased  emphasis  on  the  baro  altimeter  and  radar  altimeter  (decision  height 
(DH)  is  based  on  the  altimeters). 

•  Establishes  a  systematic  scan  for  the  runway  environment  prior  to  reaching 
DH. 

•  Continues  descent  to  DH. 

Note:  The  precision  ILS  approach  must  be  discontinued  if  the  localizer  course  becomes 
unreliable,  or  anytime  full-scale  deflection  of  the  pilot's  control  display  indicator  (CDI) 
occurs  on  final  approach  [16].  The  pilot  must  not  descend  below  localizer  minimums  if 
the  aircraft  is  more  than  one  dot  (half  scale)  below  or  two  dots  (full  scale)  above  the  glide 
slope.  If  the  glide  slope  is  recaptured  to  within  the  above  tolerance,  descent  may  be 
continued  to  DH. 

A  block  diagram  of  the  NRS  [45]  is  shown  in  Figure  1-2,  Navigation  Reference 
System  (NRS).  A  block  diagram  of  the  LSM  is  shown  in  Figure  1  -  3.  A  “walk-through” 
of  Figure  1-3  can  be  found  in  Chapter  3,  Section  3.2. 


1.9  Summary 

This  chapter  has  given  a  brief  overview  of  the  thesis  plan  to  develop  an  integrated 
GPS/INS/Baro  and  Radar  Altimeter  System  for  aircraft  precision  approaches.  Past 
research,  the  scope  of  this  project,  and  all  assumptions  were  presented.  In  Chapter  2,  the 
reference  frames  utilized  in  this  project  are  presented,  along  with  INS,  GPS,  baro 
altimeter  and  radar  altimeter  subsystems,  and  Kalman  filtering  algorithms  are  discussed. 
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Navigation  Reference  System  (NRS) 


Figure  1-2.  Navigation  Reference  System  (NRS) 
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Figure  1-3.  Landing  System  Model  (LSM)  Simulation 
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Chapter  3  introduces  and  develops  the  landing  system  model  (LSM).  Chapter  4  discusses 
the  results  and  analysis  of  the  LSM.  Chapter  5  summarizes  the  research  and  presents 
recommendations  and  conclusions  based  on  results  presented  in  Chapter  4. 


l  - 17 


II.  Background 


2.1  Introduction 

This  section  presents  basic  theory  of  a  ring  laser  gyro  (RLG),  Inertial  Navigation 
Systems  (INS),  Global  Positioning  System  (GPS),  barometric  altimeter  and  the  radar 
altimeter.  ELS  precision  approach  information  will  also  be  presented.  Fundamental 
Kalman  filter  and  extended  Kalman  filter  (EKF)  theory  will  also  be  discussed.  A  more 
rigorous  development  of  many  of  the  Kalman  filter  subjects  can  be  found  in  [40,41,42]. 

A  quick  review  of  notation  usage  in  Section  2.1.1  may  be  advantageous  before  reading 
further  in  this  chapter. 

2.1.1  Notation 

Notation  used  in  this  thesis  will  attempt  to  maintain  consistency  with  [40] 
Deterministic  and  stochastic  processes  alike  will  be  indicated  by  the  roman  typeface. 
Vectors  will  be  displayed  in  bold-faced  type,  x,  whereas  scalars  will  be  normal  type,  x. 
Matrices  will  be  displayed  in  bold-upper  case  X.  A  particular  realization  of  a  variable 
will  be  displayed  in  italics,  x. 

2.2  Ring  Laser  Gyro  (RLG)  Strapdown  INS: 

In  order  to  "navigate"  well,  one  must  find  a  solution  to  the  following  layman’s 
“navigation”  problem  [32]: 

•  Where  am  I?  (Present  position) 

•  In  which  direction  is  my  destination?  (Relative  bearing) 

•  How  far  is  it  to  my  destination?  (Distance  to  go) 

•  How  fast  and  in  what  direction  am  I  going?  (Velocity) 

A  ring  laser  gyro  strapdown  inertial  navigation  system  can  easily  find  the  solution  to  the 
above  questions. 

By  definition  [33],  an  inertial  navigator  is  a  self-contained,  dead-reckoning 
navigation  aid  using  inertial  sensors,  a  reference  direction,  and  initial  or  subsequent  fixes 
to  determine  direction,  distance,  and  speed;  single  integration  of  acceleration  provides 
speed  information  and  a  double  integration  provides  distance  information. 

From  [9],  strapdown  systems  are  characterized  by  their  lack  of  a  gimbal  support 
structure.  An  advantage  of  strapdown  vs.  the  gimbaled  is  that  a  strapdown  has  no  moving 
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platform  keeping  a  "stable  element"  level.  Eliminating  gimbals  (i.e.  reducing  moving 
parts)  translates  to  a  higher  mean  time  between  failure  (MTBF),  and  typically  lower 
production  costs.  Also,  if  a  gyro  fails  in  a  strapdown  system,  that  one  gyro  can  be 
replaced;  in  a  gimbaled  system,  the  entire  inertial  measurement  unit  (IMU)  would  have  to 
be  replaced.  A  disadvantage  strapdown  has  vs.  gimbaled  is  that  the  platform  is  physically 
“strapped-down”  to  the  dynamic  body  in  question,  thus  the  gyroscopes,  accelerometers 
and  strapdown  computer  algorithms  must  be  able  to  “keep-up”  with  whatever  harsh 
dynamic  environment  to  which  it  is  attached.  Accuracy  is  often  lower  also. 

The  strapdown  system  is  mechanized  by  mounting  three  gyros  and  three 
accelerometers  directly  to  the  vehicles  for  which  the  navigation  function  is  to  be 
provided.  An  onboard  digital  computer  keeps  track  of  the  vehicle’s  attitude  with  respect 
to  some  reference  frame  based  on  information  from  the  gyros.  The  computer  is  thus  able 
to  provide  the  coordinate  transformation  necessary  to  coordinatize  the  accelerometer 
outputs  in  a  computational  reference  frame. 

RLG  construction  typically  consists  primarily  of  an  optical  cavity,  a  laser  device, 
three  or  four  mirrors,  a  prism,  and  a  pair  of  photo  detectors  [58,54].  According  to 
[58,54,29],  the  RLG  operates  as  follows:  the  laser  gyro  detects  and  measures  angular 
rates  by  measuring  the  frequency  difference  between  two  contra-rotating  (laser)  beams 
[29].  The  two  laser  beams  circulate  in  the  "ring"  cavity  simultaneously.  Mirrors  are  used 
to  reflect  each  beam  around  an  enclosed  area.  The  resonant  frequency  of  a  contained 
laser  beam  is  a  function  of  its  optical  path  length.  Consequently,  the  two  laser  beams 
have  the  same  frequencies  when  the  gyro  is  at  rest.  If  the  cavity  is  rotating  in  an  inertial 
sense,  the  propagation  times  of  the  two  light  beams  are  different.  The  delay  manifests 
itself  in  the  form  of  a  phase  shift  between  the  two  beams,  and  the  phase  shift  is  detected 
by  a  pair  of  photo  detectors  [58,54].  The  magnitude  of  the  phase  shift  provides  a  direct 
indication  of  the  angular  rate  of  rotation  of  the  instrument  with  respect  to  inertial  space 
[58,54].  The  laser  gyro  is  an  unconventional  gyro  since  it  does  not  have  a  spinning  rotor. 
Devices  of  this  type  are  extremely  reliable  due  to  the  absence  of  moving  parts  [58,  54]. 

Specific  force  is  measured  by  accelerometers.  The  most  common  accelerometers 
to  date  have  been  devices  which  are  sophisticated  variations  of  the  simple  pendulum 
[58,9].  The  motion  of  the  internal  mass  is  related  to  the  inertially  referenced  motion  of  the 
instrument  case  by  Newton’s  second  law  of  motion  [58].  However,  to  obtain  the  correct 
measure  of  inertial  acceleration,  the  effects  of  local  gravity  must  be  removed  from  the 
measured  specific  force  [9,  52]. 

In  summary,  inertial  navigation  principles,  have  been  well  understood  for  many 
years  [9].  The  author  has  always  been  fascinated  at  the  cost  and  complexity  of  some 
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INSs,  of  how  technicians  and  engineers  continually  try  to  "tweak",  upgrade  and  massage 
its  sensor  components.  In  this  brief  INS  introduction  the  author  has  mentioned  many 
things,  one  of  them  is  that  the  INS  solves  our  navigation  problem,  and  it  does  so  on  its 
own.  A  system  self-contained,  which  does  not  need  outside  help,  cannot  be  jammed  and 
does  not  emit  radiation  that  otherwise  could  be  detected  (for  military  purposes).  The  INS 
typically  outputs  data  at  a  50  Hz  rate.  Amazingly,  INS  prices  keep  dropping  ($127K 
(gimbaled  INS)  in  1990,  now  $65K  (RLG  INS)  in  1994  -  with  double  the  reliability! 

2.3  Barometric  Altimeter: 

A  shortcoming  of  any  INS  is  the  instability  present  in  the  vertical  channel  which 
(in  the  absence  of  aiding  information)  results  in  unbounded  error  growth  in  vertical 
position  and  velocity  [3,9,25,52].  This  inherent  instability  is  controlled  by  vertical 
channel  aiding.  Such  aiding  is  frequently  accomplished  in  vertical  position  information 
provided  from  a  barometric  altimeter.  This  external  altitude  information  has  the  effect  of 
stabilizing  the  vertical  channel  [9]. 

There  are  many  ways  to  measure  the  altitude  of  an  aircraft;  probably  the  simplest 
is  with  a  barometric  altimeter  [16].  The  pressure  of  the  Earth’s  atmosphere  decreases  as 
height  above  the  earth  increases.  Barometric  altimeters  are  designed  to  output  altitude 
relative  to  the  pressure  difference.  For  example,  if  the  barometric  scale  is  referenced  to 
29.92”  Hg  (sea  level,  standard  conditions)  and  the  instrument  is  supplied  with  a  static 
pressure  of  20.58”  Hg  (pressure  at  10,000  feet,  standard  conditions),  the  altimeter  should 
output  10,000  feet.  The  pressure  difference  between  the  sea  level  and  10,000  feet  on  a 
standard  day  is  9.34  inches  of  mercury  [16],  Barometric  altimeters  are  most  inaccurate 
when  ascending  or  descending  at  rapid  rates  (especially  noted  with  fighter  aircraft)  but  are 
relatively  low  cost. 

2.4  Global  Positioning  System  (GPS) 

GPS  navigation  presents  opportunity  for  standardized  worldwide  civil  aviation 
operations  using  a  common  navigation  receiver  [23].  GPS  is  a  space-based  positioning, 
velocity  and  time  system  that  has  three  major  segments:  Space,  Control  and  User. 


2-3 


2.4.1  GPS  Space  Segment 


The  GPS  Space  Segment  is  composed  of  24  satellites  in  six  orbital  planes.  The 
satellites  operate  in  near-circular  20,200  km  (10,900  NM)  orbits  at  an  inclination  angle  of 
55  degrees  and  with  ~  12-hour  period.  The  spacing  of  satellites  in  orbit  is  arranged  so  that 
a  minimum  of  five  satellites  will  be  in  view  to  users  worldwide,  with  a  position  dilution 
of  precision  (PDOP)  of  six  or  less.  PDOP  is  a  measure  of  the  error  contributed  by  the 
geometric  relationships  of  the  GPS  satellites  as  seen  by  the  GPS  receiver  [17].  PDOP  is 
mathematically  defined  as: 

PDOP  =  (a2x+o2y+o2z)ul  (2.1) 

2  2  2 

where  Gr,o},  andaT  are  the  variances  of  the  x,  y  and  z  pseudorange  measurement 

position  errors  [17].  Each  satellite  transmits  on  two  L  band  frequencies,  LI  (1575.42 
MHz)  and  L2  (1227.6  MHz).  LI  carries  a  precise  (P)  code  and  a  coarse/acquisition  (C/A) 
code.  L2  carries  the  P  code.  A  navigation  data  message  is  superimposed  on  these  codes. 
The  same  navigation  data  message  is  carried  on  both  frequencies. 

2.4.2  GPS  Control  Segment 

The  Control  Segment  has  five  monitor  stations,  three  of  which  have  uplink 
capabilities.  The  monitor  stations  use  a  GPS  receiver  to  track  all  satellites  in  view 
passively  and  thus  accumulate  ranging  data  from  the  satellite  signals.  The  information 
from  the  monitor  stations  is  processed  at  the  Master  Control  Station  (MCS)  to  determine 
satellite  orbits  and  to  update  the  navigation  message  of  each  satellite.  This  updated 
information  is  transmitted  to  the  satellites  via  the  ground  antennas,  which  are  also  used 
for  transmitting  and  receiving  satellite  control  information. 

2.4.3  GPS  User  Segment 

The  User  segment  consists  of  an  antenna  and  receiver  processors  that  provide 
positions,  velocity  and  precise  timing  to  the  respective  user.  Computing  the  user’s 
positional  information  typically  requires  simultaneous  solution  of  the  following  four 
nonlinear  position  equations  [17]: 
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(x,  —  ux)2  +(y,  -  uy )2  +(z,  -m.)2  =  (R,  —  CB)2 
(x2  -  ux  f  +  (y2  -  My  )2  +  (Z2  -  Uz  f  =  (R2  -  CBf 

(X3  "My)2  +(y3  “My)2  +(Z3  -M^2  =(R3  -CS)2 

(X4  -  My )2  +  (y4  -  My )2  +  (z4  -uz)2=  (R4  -  CB)2 

where  the  pseudo  range,  Ri=i,2,3,4  to  each  satellite  is  defined  as 

R,  =  CAt] 

R2  =  CAt2 

r3  =  ca?3 

R4  =  CAt4 


and 


C  =  speed  of  light 

A/i=i,2,3,4  =  time  signals  transmitted  by  the  i-th  satellite 
Xj=it 2,3,4,  yi=i, 2,3,4.  zi=i,2,3,4  are  respective  i-th  satellite 
positions 

mx,  My,  uz  is  the  user  position  the  GPS  user 
equipment  is  solving  numerically  and  recursively 
CB  =  the  user  clock  bias  (user  equipment  solves) 

Normally  the  user  equipment  needs  to  acquire  and  maintain  lock  on  four  satellites 
in  order  to  compute  a  3-D  position  fix  [44]  and  the  clock  bias  CB.  The  GPS  pseudorange 
between  the  user  and  each  satellite  is  computed  based  on  knowledge  of  time  (the  master 
GPS  clock)  and  the  unique  signal  format  which  is  broadcast  by  each  satellite.  Once  the 
four  pseudo-ranges  are  known,  a  recursive  algorithm  is  solved  to  compute  the  user’s 
position  [44],  See  [17]  for  further  references. 

2.5  Radar  Altimeter 

A  radar  altimeter  provides  measurement  of  absolute  clearance  over  all  types  of 
terrain  [28].  System  operation  is  based  on  the  precise  measurement  of  the  time  required 
for  an  electromagnetic  energy  pulse  to  travel  from  the  aircraft  to  the  terrain  below  and  to 
return.  An  elementary  block  diagram  (transmitter,  receiver,  range  computer,  height 
indicator)  is  shown  in  Figure  2  -  1 . 
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Figure  2-1.  Elementary  Radar  Altimeter  Block  Diagram 


Radar  altimeters  are  normally  all-weather  devices.  Performance  specifications  (3- 
a)  are  typically  ±  [3ft  +  3%  of  altitude  range],  with  ±  30°  pitch  and  ±  45°  roll 
maneuverability  at  above  ground  level  (AGL)  heights,  which  typically  vary  from  0  feet  to 
10,000  feet. 

2.6  Instrument  Landing  System  (ILS)  Precision  Approach 

An  instrument  approach,  by  definition  [33],  is  the  process  of  making  an  approach 
to  a  landing  by  the  use  of  navigation  instruments  without  dependence  upon  direct  visual 
reference  to  the  terrain.  The  Instrument  Landing  System  (ILS)  is  designed  to  provide  an 
approach  path  for  exact  alignment  and  descent  of  an  aircraft  on  final  approach  to  a 
runway  [16].  The  ground  equipment  consists  of  two  highly  directional  transmitting 
systems,  and  along  the  approach,  three  (or  fewer)  marker  beacons.  The  directional 
transmitters  are  known  as  the  localizer  and  glide  slope  transmitters.  See  Figure  2-2. 

The  localizer  transmitter,  operating  on  one  of  the  40  ILS  channels  within  the 
frequency  range  of  108. 10  MHz  to  1 1 1.95  MHz,  emits  signals  which  provide  the  pilot 
with  course  guidance  to  the  runway  centerline  in  the  horizontal  plane.  The  localizer 
signal  is  usable  and  accurate  to  a  range  of  18  nautical  miles  (NM)  from  the  localizer 
antenna  unless  otherwise  depicted  on  the  Instrument  approach  procedure  (LAP)  [16].  See 
Figure  2-2. 

The  ultra  high  frequency  glide  slope  transmitter,  operating  on  one  of  the  40  ILS 
channels  within  the  frequency  range  329.15  MHz  to  335.00  MHz,  radiates  its  signals 
primarily  in  the  direction  of  the  localizer  front  course,  i.e.,  so  as  to  measure  angular 
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vertical  displacement  from  the  desired  glide  path,  as  seen  from  the  side.  The  glide  slope 
signal  is  usable  to  a  distance  of  10  NM  unless  otherwise  depicted  on  the  Instrument 
approach  procedure  (IAP)  [16]. 

A  marker  beacon  light  and  (or)  aural  tone  may  be  included  in  the  cockpit  display 
to  indicate  aircraft  position  along  the  localizer.  The  marker  beacons  are  identified  by 
continuous  dashes  for  the  outer  marker,  alternating  dashes  and  dots  for  the  middle 
marker,  and  continuous  dots  for  the  inner  marker.  See  Figure  2-2.  Precision  ILS 
Approaches  follow  the  2-a  decision  height  accuracy  shown  in  Table  2-1,  2-cr  Precision 
ILS  Approach  Criteria  at  Decision  Height. 


ILS  Precision  Category 

Horizontal  Accuracy 

Vertical  Accuracy 

1 

±56.1  ft 

±  13.5  ft 

2 

±  17.1  ft 

±5.6  ft 

3 

±  13.5  ft 

‘±  2.0  ft 

Table  2-1.  2-g  Precision  ILS  Approach  Criteria  at  Decision  Height 


2.7  Reference  Frames 

A  navigation  "solution"  has  significance  only  if  the  corresponding  frame  in  which  the 
solution  is  expressed  is  clearly  understood  [58].  While  the  preceding  statement  may  seem 
obvious,  it  cannot  be  overemphasized.  Consider  that  a  typical  INS  "owner’s  manual" 
defines  earth  frame,  true  frame,  computer  frame,  platform  frame,  sensor  frame, 
accelerometer  frame  and  the  body  frame  [58,38].  From  a  student’s  perspective  this  may 
at  first  be  overwhelming,  but  to  make  matters  worse,  another  INS  vendor  may  well  define 
every  frame  mentioned  above,  such  as  "earth  frame"  in  an  entirely  different  manner! 
Therefore,  the  frames  used  in  this  project,  along  with  coordinate  transformations,  will 
briefly  be  discussed  (all  reference  frame  figures  with  permission  from  [7]). 
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Figure  2  -  2.  Standard  ILS  Characteristics  and  Terminology. 
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2.7.1  Inertial  Frame  {x,  y\z) 

An  inertial  frame  is  an  orthogonal,  right-handed  coordinate  system;  its  origin  is 
coincident  with  the  earth’s  center-of-mass  and  the  frame  is  oriented  as  follows.  The  x‘,  z 
plane  lies  in  the  earth’s  equatorial  plane  and  does  not  rotate  with  respect  to  the  fixed  stars. 
The  y1  axis  projects  from  the  earth’s  center-of-mass  directly  through  the  North  pole.  The 
inertial  frame  is  depicted  by  the  [x1,  y1,  z]  frame  shown  in  Figure  2-3. 

2.7.2  Earth  Frame  (xe,  ye,  ze) 

The  earth  frame  or  “earth-centered-earth-fixed”  (ECEF)  frame  is  an  orthogonal, 
right-hand  coordinate  system;  its  origin  is  coincident  with  the  earth’s  center-of-mass, 
with  the  xe,  z  plane  located  in  the  earth’s  equatorial  plane.  The  z  axis  is  aligned  with  the 
Greenwich  meridian  and  rotates  at  exactly  the  earth  rate,  Q,  about  the  ye  axis,  which 
projects  from  the  earth’s  center-of-mass  directly  through  the  North  pole.  The  Earth  frame 
is  depicted  as  [xe,  ye,  ze]  in  Figure  2-4. 

2.7.3  Geographic  Frame  (x8,  y8,  z8)  =  (E,  N,  U) 

The  geographic  frame  or  “local-level”  frame  is  an  orthogonal,  right-handed 
coordinate  system;  its  origin  is  at  the  location  of  the  INS  (or  the  user),  and  its  axes  are 
aligned  with  the  East,  North  and  Up  directions  [E,  N,  U].  The  geographic  frame  remains 
perpendicular  to  the  earth’s  surface  with  respect  to  the  earth’s  gravity  field  as  the  user 
moves  over  the  Earth.  The  geographic  frame  is  depicted  as  either  [xg,  yg,  zg]  or  [E,  N,  U] 
in  Figure  2-5. 

2.7.4  Navigation  Frame  (xn,  yn,  zn) 

The  navigation  frame  or  “local-level-wander-azimuth”  frame  is  an  orthogonal, 
right-hand  coordinate  system;  its  origin  is  at  the  location  of  the  INS  (or  the  user).  This 
frame  coincides  with  the  geographic  frame  when  the  wander  angle,  a  equals  0°.  The 
wander  angle  is  a  computed  angle  between  a  “scribe  mark”  on  a  wander  azimuth  angle 
platform  and  North.  For  gimbaled  systems,  the  platform  is  purposely  not  commanded  to 
seek  North,  due  to  the  high  platform  angular  rates  that  this  would  require  in  polar  regions, 
with  resulting  performance  degradation  [3,58].  The  navigation  frame  is  denoted  as  [xn, 
yn,  zn]  shown  in  Figure  2-6. 
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INERTIAL  FRAME  (x*.  y',  z‘) 


Figure  2-3.  Inertial  Frame 


EARTH  FRAME  tx*,  y*,  z*) 


Figure  2-4.  Earth  Frame 


$ 


GEOGRAPHIC  FRAME  <x«,  y»,  z®)  -  |E,  N,  U) 


PLANE 


Figure  2-5.  Geographic  Frame 


LOCAL  LEVEL,  WANDER  AZIMUTH,  NAVIGATION  FRAME  (xn.  y",  z") 


PLANE 


Figure  2-6.  Navigation  Frame 


2.7.5  Body  Frame  (xb,  yb,  zb) 

The  body  frame  is  an  orthogonal,  right-hand  frame;  its  origin  is  at  the  vehicle  (i.e., 
aircraft)  center-of-mass.  Its  axes  are  the  vehicle’s  roll,  pitch,  and  yaw  axes  [<j>,  6,  vj/].  The 
xb  axis  points  in  the  forward  direction,  along  the  roll  axis;  the  yb  axis  points  to  the  right 
(starboard  side)  of  the  aircraft,  perpendicular  to  the  roll  axis,  but  along  the  pitch  axis;  and 
the  zb  axis  is  positive  out  the  underside  of  the  aircraft.  The  body  frame  is  denoted  as  [xb, 
yb,  zb]  and  is  shown  in  Figure  2-7. 

2.8  Reference  Frame  Transformations 

The  RLG  INS  modeled  in  this  thesis  uses  the  navigation  frame  or  “local-level- 
wander-azimuth”  frame.  It  is  often  necessary  to  express  vectors  such  as  position,  attitude, 
velocity  or  acceleration  in  terms  of  several  different  reference  frames.  As  an  example,  the 
INS  modeled  in  this  thesis  also  outputs  position  error  in  terms  of  an  error-angle  vector, 
[50x,  50y,  80z,  Sh]T,  where  80x  is  the  error  angle  about  the  local  level  xg  (or  E)  axis,  S0y  is 
the  error  angle  about  the  local  level  yg  (or  N)  axis,  80z  is  the  error  about  the  local  level  zg 
(or  U)  axis,  and  8h  is  the  altitude  error  [3,38].  Even  though  this  definition  is  clear,  if  the 
error-angle  vector  is  to  have  physical  meaning,  it  must  be  transformed  into  a  vector  in 
navigation  error  terminology,  [&)),  8X,  8a,  Sh]T,  where  8$  is  the  error  in  latitude,  8X  is 
longitude  error,  8a  is  alpha  angle  error  and  8h  is  again  the  altitude  error. 


BODY  FRAME  (xb  yb,  zb 
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A  transformation  matrix,  C^ro? Angle Error ,  permits  compact  transformation  of  the 


error-angle  vector  into  an  equivalent  expression  in  navigation  error  space.  The 
transformation  matrix,  C!^^fa^leError  is  shown  below  in  Equation  (2-1)  [3,58]: 


-  cos  a 

^Navigation  Error  =  sin  «  SeC  <l> 
'-Error  Angle  -  sin  a  tan  <|> 

0 


sin  a 
cos  a  sec  <j) 
-cos  a  tan  (j) 
0 


0 

0 

1 

0 


Other  transformations  are  as  follows  (from  [7]): 
2.8.1  Inertial  Frame  to  Earth  Frame,  Cf 


Cf  = 


cos Qt  0  -sin  Qt 
0  1  0 
sin  Clt  0  cos  Q.t 


2.8.2  Earth  Frame  to  Geographic  Frame,  Cf 


0 

0 

0 

1 


(2.2) 


(2.3) 


'1 

0 

0 

cos  A, 

0 

-sin  A, 

0 

COS  <|) 

-sin  <J) 

0 

1 

0 

0 

sin  (j) 

COS<J) 

sin  A, 

0 

cos  A, 

cos  A,  0  -sin  A, 

-  sin  <j>  sin  A,  cos  <)>  —  sin  c}>  cos  A. 

cos  (|)  sin  A.  sin  <|)  cos  (j)  cos  A, 


where 


(2.4) 


A,  =  longitude 
(f>  =  latitude 
a  =  alpha  angle 
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2.8.3  Earth  Frame  to  Navigation  Frame, C” 


cos  a 

sin  a 

0 

1 

0 

0 

cos  A. 

0 

-sin  A, 

cn  = 

-sin  a 

cos  a 

0 

0 

COS(|) 

-sintj) 

0 

1 

0 

0 

0 

1 

0 

sintj) 

COStj) 

sin  A. 

0 

cos  A. 

cos  a  cos  A, -sin  a  sin  <j)  sin  A.  sinacostj)  -  cos  a  sin  A, -sin  a  sin  <{)  cos  A, 

-sinacosA -  cos  a  sin  (j)  sin  A  cosacostj)  sin  a  sin  A  -  cos  a  sin  (j)  cos  A 

cosp  sin  A  sin(j)  cos  p  cos  A 


where 


A  =  longitude 
<|)  =  latitude 
a  =  alpha  angle 

2.8.4  Geographic  Frame  to  Navigation  Frame,  C" 

6 


where 


C" 


cos  a  sin  a  0 
-sin  a  cos  a  0 
0  0  1 


a  =  alpha  angle 


2.8.5  Geographic  Frame  to  Body  Frame,  Cbp 

~  o 


'1 

0 

0 

COS0 

0 

-sin0 

cosy 

siny 

O' 

'0 

1 

0  ' 

Cb  = 
C 8 

0 

cosp 

sinp 

0 

1 

0 

-siny 

cosy 

0 

1 

0 

0 

0 

-sinp 

cosp 

sin0 

0 

COS0 

0 

0 

1 

0 

0 

-1 

cos  0  sin  \|/  cos  0  cost)/  sin0 

sin p sin 0 sin  \|/  + cos p cosy  sin p sin 0 cos y  -  cos p sin y  -sinpcos0 
cos  p  sin  0  sin  y- sin  p  cosy  cospsin0cosy +  sinpsiny  -cospcos© 


where 

p  =  roll 
0  =  pitch 

y  =  geographic  heading 


(2.6) 


(2.7) 
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L. 

2.8.6  Navigation  Frame  to  Body  Frame,  Cn 


"1 

0 

0 

COS0 

0 

-sin0 

COSVj/p 

sin  \\f  P 

O' 

'0 

1 

O' 

cb  = 

0 

cosp 

sinp 

0 

1 

0 

-  sin  t]/  P 

cost]/  p 

0 

1 

0 

0 

0 

-sinp 

cosp 

sin0 

0 

cos0 

0 

0 

1_ 

_0 

0 

-1 

cosSsinij/p  cos0cosv)/p 

sin  p  sin  0  sin  \|/ P  +cospcos  \\t  P  sinp  sin  0  cos  tj/p  -cosp  sin  t|//> 
cosp  sin  0  sin  \|/P  -sinpcost|/P  cosp  sin  0  cos  tj/P  +  sinpsint|/P 


sin0 

-sinpcos© 

-cospcos0 


-,(2.8) 


where 


p  =  roll 
0  =  pitch 

\|/  P  =  platform  heading 


2.9  Kalman  Filter  Theory 
2.9.1  What  is  a  Kalman  Filter? 

A  Kalman  filter  is  simply  an  optimal  recursive  data  processing  algorithm  [40]  that 
can  be  shown  to  be  optimal  by  essentially  any  standard,  given  the  appropriateness  of 
several  underlying  assumptions.  These  assumptions  are  that  the  system  in  question  can 
be  adequately  modeled  as  linear  with  white,  Gaussian  system  and  measurement  noises. 

One  aspect  of  the  word  "optimal"  is  that  the  Kalman  filter  can  incorporate  all 
information  (measurements)  provided  to  it  [40].  It  processes  all  available  measurements 
regardless  of  their  precision  to  "estimate"  the  current  value  of  the  variables  of  interest 
with  use  of  (from  [40]): 

•  Knowledge  of  the  system  and  measurement  device  dynamics 

•  The  statistical  description  of  the  system  noises,  measurement  errors  and 
uncertainty  in  the  dynamics  models. 

•  Any  available  information  about  inertial  conditions  of  the  variables  of 
interest. 

For  example,  to  determine  the  velocity  of  an  aircraft,  one  could  use  a  Doppler 
radar,  or  the  velocity  indications  from  an  inertial  navigation  system,  or  the  pitot  and  static 
pressure  and  relative  wind  information  in  the  air  data  system.  Rather  than  ignore  any  of 
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these  outputs,  a  Kalman  filter  could  be  built  to  combine  all  this  data  and  knowledge  of  the 
various  systems  dynamics  to  generate  an  overall  best  estimate  of  velocity.  Another  way  a 
Kalman  filter  is  optimal  is  that  it  obtains  the  best  estimate  of  desired  quantities  from  data 
provided  by  a  noisy  environment.  Here  the  word  “optimal”  means  that  the  Kalman  filter 
minimizes  errors  in  essentially  all  respects,  and  it  does  so  recursively.  The  word 
recursive  means  that,  unlike  certain  data  processing  concepts,  the  Kalman  filter  doesn't 
require  all  previous  data  to  be  kept  in  storage  and  reprocessed  every  time  new 
measurements  are  taken. 

To  “see”  how  a  Kalman  filter  works,  a  simple  example  taken  directly  from  [39] 
will  be  presented.  It  is  included  it  here  because  it  helped  the  author  understand  the 
concept  of  a  Kalman  Filter  in  his  AFIT  studies. 

2.9.2  Kalman  Filter  Example 

Suppose  that  you  are  lost  at  sea  during  the  night  and  have  no  idea  at  all  of  your 
location.  So  you  take  a  star  sighting  to  establish  your  position  (for  the  sake  of  simplicity, 
consider  a  one-dimensional  location).  At  some  time  r,  you  determine  your  location  to  be 
Z\ .  However,  because  of  inherent  measuring  device  inaccuracies,  human  error,  and  the 

like,  the  result  of  your  measurement  is  somewhat  uncertain.  Say  you  decide  that  the 
precision  is  such  that  the  standard  deviation  (one-sigma  value)  involved  is  a  (or 


Figure  2-8.  Conditional  Density  of  Position  Based  on  Measured  Value  zx 


2-16 


2 

equivalently,  the  variance,  or  second  order  central  statistic,  is  aZ| ).  Thus,  you  can 
establish  the  conditional  probability  of  x (tt ) ,  your  position  at  time  tx ,  conditioned 
on  the  observed  value  of  the  measurement  being  Z\ ,  as  depicted  in  Figure  2-8.  This  is  a 
plot  of  fx(t  ^  (x\z,\  )as  a  function  of  the  location  x:  it  tells  you  the  probability  of  being 

in  any  one  location,  based  upon  the  measurement  you  took.  Note  that  is  a  direct 
measure  of  the  uncertainty:  the  larger  aZ|  is,  the  broader  the  probability  peak  is,  spreading 

the  probability  "weight"  over  a  larger  range  of  x  values.  For  a  Gaussian  density,  68.3%  of 
the  probability  "weight"  is  contained  within  the  band  a  units  to  each  side  of  the  mean,  the 
shaded  portion  in  Figure  2-8. 

Based  this  conditional  probability  density,  the  best  estimate  of  your  position  is 

*(r,)  =  zi  (2.9) 

and  the  variance  of  the  error  in  the  estimate  is 


=  of,  (2.10) 

Note  that  x  is  both  the  mode  (value  that  locates  the  peak)  and  the  median  (value  with  1/2 
of  the  probability  weight  to  each  side),  as  well  as  the  mean  (center-of-mass). 

Now  say  a  trained  navigator  friend  takes  an  independent  fix  right  after  you  do,  at 
time  t2  =  ty  (so  that  the  true  position  has  not  changed  at  all),  and  obtains  a  measurement 
z2  with  a  variance  aZ2 .  Because  he  has  a  higher  skill,  assume  the  variance  in  his 

measurement  to  be  somewhat  smaller  than  in  yours.  Figure  2-9  presents  the  conditional 
density  of  your  position  at  time  t2,  based  only  on  the  measured  value  z,2.  Note  the 
narrower  peak  due  to  smaller  variance,  indicating  that  you  are  rather  certain  of  your 
position  based  on  his  measurement. 

At  this  point,  you  have  two  measurements  available  for  estimating  your  position. 

The  question  is,  how  do  you  combine  these  data?  It  can  be  shown  that,  based  on  the 
assumptions  made,  the  conditional  density  of  your  position  at  t2  =  tx ,  x(t2  )given  both  Z\ 
and  z2,  is  a  Gaussian  density  with  mean  \i  and  variance  a2  as  indicated  in  Figure  2-10 
with 


P-  = 


+  o: 


Z\  + 


at  +  at 


*2 


(2.11) 
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Figure  2-9.  Conditional  Density  of  Position  Based  on  Measurement  z2  Alone 


Figure  2-10.  Conditional  density  of  position  based  on  data  zj  and  z2 
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1 


(2.12) 


Note  that,  from  (2. 12),  o  is  less  than  either  a  or  aZi ,  which  is  to  say  that  the  uncertainty 

in  your  estimate  of  position  has  been  decreased  by  combining  the  two  pieces  of 
information. 

Given  this  density,  the  best  estimate  is 

x(f2)  =  M-  (2-13) 


with  an  associated  error  variance  a  .  It  is  the  mode,  the  median  and  the  mean  (or,  since  it 
is  the  mean  of  a  conditional  density,  it  is  also  termed  the  conditional  mean).  Furthermore, 
it  is  also  the  maximum  likelihood  estimate,  the  weighted  least  squares  estimate,  and  the 
linear  unbiased  estimate  whose  variance  is  less  than  that  of  any  other  linear  unbiased 
estimate.  In  other  words,  it  is  the  "best"  you  can  do  according  to  just  about  any 
reasonable  criterion. 

After  some  study,  the  form  of  g  given  in  (2.1 1)  makes  good  sense.  If  GZ|  were 
equal  to  ctZj  ,  which  is  to  say  you  think  the  measurements  are  of  equal  precision,  the 

equation  says  the  optimal  estimate  of  position  is  simply  the  average  of  the  two 
measurements,  as  would  be  expected.  On  the  other  hand,  if  GZ|  were  larger  than  gZj  , 

which  is  to  say  that  the  uncertainty  involved  in  the  measurement  Z\  is  greater  than  that  of 
Z2,  then  the  equation  dictates  "weighting"  z2  more  heavily  than  Z\.  Finally,  the  variances 
of  the  estimate  is  less  than  GZ|  even  if  oZjis  very  large:  even  poor  quality  data  provides 

some  information,  and  should  thus  increase  the  precision  of  the  filter  output. 

The  equation  for  x(t2  )can  be  rewritten  as 


x(t2)  = 


r  2 

r  ?  i 

<37 

z2 

Zi  + 

2  2 

2  ,  2 

o:  +  o: 

L  zt  z2  J 

L  Gl  +  °z2  \ 

Z2 


(2.14) 


or 


x(t2)  =  z{  + 


[z2  -  Zi  ] 


(2.15) 


or,  in  final  form  that  is  actually  used  Kalman  filter  implementations  (noting  that 

*(?i)  =  Zi), 
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x(t2)  =  x(tl)  +  K(t2)[z2-x(t01 


(2.16) 


where 


K{t2) 


(2.17) 


These  equations  say  the  optimal  estimate  at  t2,  x(t2  ),  is  equal  to  the  best  prediction  of  its 
value  before  z2  is  taken,  x  ( f , ) ,  plus  a  correction  term  of  an  optimal  weighting  value 
times  the  residual  difference  between  z2  and  the  best  prediction  of  its  value  before  it  is 
actually  taken,  x  ( r, ) .  It  is  worthwhile  to  understand  this  "predictor-corrector"  structure 
of  the  filter.  Based  on  all  previous  information,  a  prediction  of  the  value  that  the  desired 
variables  and  measurement  will  have  at  the  next  measurement  time  is  made.  Then,  when 
the  next  measurement  is  taken,  the  difference  between  it  and  its  predicted  value  is  used  to 
"correct"  the  prediction  of  the  desired  variables. 

Using  the  K(t2)  in  Equation  (2.17),  the  variance  equation  given  by  (2.12)  can  be 
rewritten  as 


=  o2x{h)~  K{t2Wx(tx)  (2.18) 

Note  that  the  values  of  x(t2  )  and  (t2)  embody  all  the  information  in 

fx  ( ,2 1  z ( ) ,  z  ( t2 )  ( •*  I  zi ,  z2 )  •  Stated  differently,  by  propagating  these  two  variables,  the 

conditional  density  of  your  position  at  time  t2,  given  Z\  and  z2,  is  completely  specified. 

Thus  we  have  solved  the  (static)  estimation  problem  example  from  [40].  This  will 
be  of  vital  importance  to  the  practicality  of  filter  implementation.  The  filter  is  plain  and 
simple,  “just  a  computer  program  in  a  central  processor”  [40].  If  the  reader  needs  a 
further  example  detailing  dynamics  and  propagations,  see  [40]. 

2.9.3  Linear  Kalman  Filter: 

Whenever  possible,  a  system  will  be  modeled  as  a  set  of  linear  differential 
equations  of  the  form  [40]: 

x(t)  =  F  (t)x(O  +  B(r)u(0  +  G(t)w(t)  (2.19) 

where: 

x  =  "state”  vector  (n-dimensional) 
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F  =  homogenous  state  dynamics  matrix  (n  x  n) 

B  =  control  input  matrix  (n  x  r) 
u  =  deterministic  control  input  vector  (  r-dimensional) 

G  =  driving  noise  input  matrix  (n  x  s) 
w  =  white  Gaussian  driving  noise  vector  (s-dimensional) 

Because  the  deterministic  control  term  is  zero  in  this  research,  it  will  be 

ignored  hereafter.  The  expected  value  (i.e.  mean),  of  the  white  Gaussian  driving  noise 
vector,  w(t)  is: 


£{w(0}  =  0 


(2.20) 


and  the  noise  strength  is  Q(t): 

£{w(0wr(f+T)}  =  Q(f)8(T)  (2.21) 

where  §(■)  is  the  Dirac  delta  function. 

While  Equation  (2.19)  is  written  in  terms  of  "whole"  value  state  variables,  the 
models  used  in  the  thesis  are  those  of  error  states.  This  choice  of  state  variable  results  in 
simpler  dynamics  equations  [9],  and  (2.19)  may  be  rewritten  as  [40]: 

5 x(t)  =  F(f)8x(f)  +  B(f)u(f)  +  G(t)w(0  (2.22) 

where  x(f)  has  been  replaced  by  the  error  state  vector  8x(f),  and  all  other  quantities  retain 
their  previous  definitions.  The  topic  of  error  states  is  explored  more  fully  in  the  section 
on  extended  Kalman  filters. 

As  previously  stated,  the  Kalman  filter  incorporates  sampled-data  measurement 
information  from  external  measuring  devices.  Irrespective  of  the  type  of  measuring 
device,  the  equation  which  is  used  to  describe  linear  measurements  is  of  the  form  [3]: 


ZG;)  =  H(tl)x(tl)  +  v(ff)  (2.23) 

or,  in  the  case  of  error-state  models: 

5z  (ff )  =  H(f,-  )S  x(f,- )  +  v(ff )  (2.24) 

where,  in  both  cases  above,  H  is  the  observation  matrix,  and  v  is  a  discrete-time  zero- 
mean  white  Gaussian  measurement  noise  vector  with  covariance  given  by  [40]: 
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(2.25) 


R(tf)  for  ti  = 
0  for  f;  *  tj 


The  Kalman  filter  "propagates"  the  error  state  and  its  covariance  from  the  instant 
in  time  immediately  following  the  most  recent  measurement  update,  t+,  to  the  instant  in 
time  immediately  preceding  the  next  measurement  update,  h+ 1’  by  numerical  integration 
of  the  following  equations  [40]: 

i(f  /  /,•)  =  ¥(t)x(t  /  /,.)  (2.26) 

P (//*,.)  =  F(r)P(?/rt)  +  P(r/rt)Fr(0  +  G(0Q(0Gr(0  (2.27) 

The  notation  for  x(t  /  q  )and  associated  error  covariance  P(t  /  tt )  indicate  the  best 
estimate  of  x  and  P  at  time  t,  based  on  measurements  through  time  q .  Initial  conditions 
are  given  as 


II 

(2.28) 

p(ff  ih)  =  p{tt) 

(2.29) 

as  provided  by  the  measurement  update  cycle  at  time  tt .  The  variables  tt  and  ti+l 

indicate  the  initial  and  final  times  for  each  integration  period,  respectively. 

After  propagation,  x(ti+l  /  ti)  =  x(r,-+1)  and  P(f/+1  /  tt)  =  P(f,-+1)are 

"updated"  (meaning  that  state  estimates  are  revised,  based  on  new  measurement 
information).  The  pivotal  element  in  the  update  equations  for  sample  time  q  shown 
below  is  the  time-varying  Kalman  filter  gain  K(tt ).  The  K(t( )  matrix  assigns  "weights" 
to  the  "measurement  residual”  (the  residual  consists  of  the  difference  between  the  actual 
measurement  and  the  filter's  prediction  of  the  measurement)  based  on  known 
measurement  noise  statistics  and  filter-computed  state  error  covariance  from  the  previous 
time  step.  This  process  is  designed  to  improve  the  estimate  of  each  element  of  the  state 
vector.  The  update  equations  are  [40]: 

K  (tt )  =  P(ff  )Hr  {f  ){H(r,.  )P(rf  )Hr  (f )  +  R(q  )}_1  (2.30) 

x(t? )  =  x(tf )  +  K (/,.  ){zf  -  H )x(ti )]}  (2.31) 

P(?I+)  =  P(^)-K(?i)H(fi)P(fr)  (2.32) 
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Although  the  algorithm  shown  above  is  generally  applicable  to  any  problem 
which  lends  itself  to  a  Kalman  filtering  solution,  it  is  not  necessarily  the  algorithm  which 
is  used  in  practice.  It  is  often  advantageous  to  use  a  form  of  the  algorithm  known  as  the 
U-D  factorization  form  [40].  In  the  U-D  algorithm,  the  filter  covariance  matrix  is  not 
propagated  as  a  single  square  array.  The  U  and  D  matrices  below  representing  the  pre- 
and  post-measurement  filter  covariances,  respectively,  are  explicitly  computed  instead 
[40]: 


p  an  =  uarwouur) 

(2.33) 

P(*+)  =  U(r/f)D(t,+  )U(tI+) 

(2.34) 

where  the  U  matrices  are  upper  triangular  and  unitary  (and  thus  contain  ones  along  the 
main  diagonal),  and  the  D  matrices  are  simply  diagonal  [40].  This  form  offers  several 
advantages  including  numerical  stability,  improved  precision,  and  guaranteed  non¬ 
negativity  of  the  computed  covariance's  eigenvalues  [40].  It  is  the  U-D  form  of  the 
Kalman  filter  algorithm  which  is  implemented  in  the  MSOFE  software  [48]  that  is  used 
in  this  research. 

2.9.4  Linearized  and  Extended  Kalman  Filtering 

Unfortunately,  not  all  problems  are  adequately  described  with  linear  systems 
driven  by  white  Gaussian  noise.  In  many  cases,  the  most  appropriate  model  is  nonlinear. 

The  navigation  problem  at  hand  falls  squarely  into  the  nonlinear  category.  Fortunately,  a 
method  exists  whereby  a  nonlinear  system  may  be  treated  in  much  the  same  manner  as  a 
linear  one  for  a  particular  class  of  problems.  Suppose  that  the  nonlinear  system  may  be 
described  by  [41]: 

x(r)  =  f  [x(r),  u(r),  t]  +  G(r)w(r)  (2.35) 

In  this  case,  the  state  dynamics  vector,  f  [., ., .],  is  a  nonlinear  function  of  the  state 
vector  x(-),  time  t,  and  the  control  input  (assumed  to  be  zero  in  this  research).  The  white 
Gaussian  noise  is  defined  exactly  as  in  (2.20)  and  (2.21),  and  it  still  enters  the  dynamics 
model  linearly.  In  addition,  the  measurement  equation  may  also  be  a  nonlinear  function 
of  the  state  vector  and  time  [41]: 


z(rt)  =  htxlr,-),*,]  +  v(rf) 


(2.36) 
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The  noise  vector  v  is  again  zero-mean,  white  and  Gaussian,  entering  the  measurement 
equation  linearly,  and  its  covariance  is  described  by  (2.25). 

Recalling  that  a  system  must  be  linear  in  order  to  satisfy  the  assumptions  upon 
which  a  conventional  Kalman  filter  is  based,  the  nonlinear  equations  (2.35)  and  (2.36) 
must  be  linearized.  The  following  approach  is  summarized  from  [41]: 

1.  Assume  that  a  nominal  state  trajectory,  xn  (t),  may  be  generated  which  satisfies 

xn(^o)  =  x„0  (2.37) 

and 

xn(r)  =  f[x„(0,u(f),f]  (2.38) 

where  f  [., ., .]  is  specified  in  (2.35),  and  u(t)=0. 

2.  The  “nominal”  measurements  which  accompany  the  nominal  trajectory  are: 


z  „(*/)  =  h[x„  (*,■),*,■] 


(2.39) 


3.  The  “perturbation”  of  the  state  derivative  is  obtained  by  subtracting  the  nominal 
trajectory  from  the  original  nonlinear  equation: 

[x(0-*„(0]  =  f[x(/),u(t),l]  -  f[x„(t),u(t),t]  +  G(t)w(t)  (2.40) 


4.  The  equation  above  may  be  approximated  to  first  order  by  a  Taylor  series  expansion: 

Sx(r)  =  F[t;xn(0]5x(0  +  G(r)w(r)  (2.41) 


where  8x(  )  represents  a  first-order  approximation  of  the  process  [x(-)  -  xn  (•)],  and 
F[r;  x„  (0]  is  a  matrix  of  partial  derivatives  of  f  with  respect  to  its  first  argument, 
evaluated  along  the  nominal  trajectory  [41]: 


F[f;x„(f)] 


3f[x,r] 


3x 


x  =  x„(r) 


(2.42) 


5.  The  perturbation  measurement  equation  is  derived  in  like  fashion  and  is  expressed  as 
[41]: 


5z(r,)  =  H[rI;x„(rI)]8x(r)+v(rI) 


(2.43) 


2-24 


where 


H[f,-;x„  (/;)]  = 


3h[x,rt] 

dx 


x  =  xn  (tt ) 


(2.44) 


With  the  “error-state”  model  in  hand,  it  is  possible  to  return  to  the  linear  filtering  theory. 

An  estimate  of  the  whole- valued  quantities  of  interest  is  obtained  from  [41]: 

x  =  xn(t)  +  8x(t)  (2.45) 


The  expression  above  for  the  linearized  Kalman  filter  is  useful,  provided  that  the 
linearization  assumption  is  not  violated.  However,  if  the  nominal  and  “true”  trajectories 
differ  by  too  large  an  amount,  unacceptable  errors  may  result  [41].  It  is  for  this  reason 
that  extended  Kalman  filtering  is  useful  in  many  cases  for  which  perturbation  techniques 
alone  do  not  suffice.  Extended  Kalman  filtering  allows  for  relinearizing  about  newly 
declared  nominals  at  each  sample  time,  to  enhance  the  adequacy  of  the  linearization 
process,  and  thus  of  the  resulting  filter  performance  as  well  [41]. 

The  extended  Kalman  filter  equations  are  summarized  below.  The  reader  is 
referred  to  [41]  for  details  regarding  their  derivation.  The  assumed  measurement  model 
equation  for  an  extended  Kalman  filter  development  is  given  by  Equation  (2.36),  where 
v(  )  is  once  again  zero-mean,  white  and  Gaussian,  with  covariance  given  by  (2.25). 
Measurements  are  incorporated  into  the  extended  Kalman  filter  via  the  following  set  of 
equations  [41]: 

K(f,- ) = par  )Hr  [tt  ;iar  )]{H[rf  ;iar  )]par  )Hr  a,  ;&ar )] + Ra,  )}_I  (2.46) 

*a,+  )  =  iafj  +  KU^z,.  -h[x(rf );;,.]}  (2.47) 


par)  =  p(tr)-K(tl)H[t,.;x(tr)]P(?r) 


where 


Ha,.;iar)] 


dh[x,f,-] 

dx 


X  =  x(ti  ) 


(2.48) 


(2.49) 


The  state  estimate  and  covariance  are  propagated  from  tt  to  ti+l  by  integrating  the 
following  equations  [41]: 


2-25 


x(t  / 1,)  =  f[x(r  /  r;),u(0,r] 

(2.50) 

p  (t/ti) 

=  F[r; x(t  /  t,)]P(t  /ti)  +  P (t  /  t!  )Fr[r;x(r  / t,-)]  +  G(OQ(OGr 

(0(2.51) 

where 

F [r,x(t/t:)]=  _  . 

(2.52) 

dx  |x  =  x(r  /  tt) 

and  the  initial  conditions  are: 

x(ti  /  /,•)  =  x(f* ) 

(2.53) 

P  (ti/ti)  =  P(t,+  ) 

(2.54) 

The  equations  shown  above  for  the  extended  Kalman  filter  are  programmed  into 
the  MSOFE  shell  [48]  for  the  problem  defined  by  this  thesis.  It  is  the  fact  that  the 
extended  Kalman  filter  is  relinearized  about  each  successive  estimate  of  the  state 
x(t)  which  “enhances  the  validity  of  the  assumption  that  deviations  from  the  reference 
(nominal)  trajectory  are  small  enough  to  allow  linear  perturbation  techniques  to  be 
employed”  [41]. 

2.10  Summary 

This  chapter  has  presented  the  basic  theory  of  RLG  INS,  GPS,  and  barometric 
altimeter  and  radar  altimeters.  ILS  precision  approach  has  also  been  defined.  Reference 
frames  and  coordinate  transformations  used  in  this  thesis  have  also  been  defined,  as  well 
as  describing  the  Kalman  filter  by  an  example  from  [40]  and  noting  linear,  linearized  and 
extended  Kalman  filter  fundamentals.  Chapter  3  will  describe  the  design  methodology 
and  error  models  of  the  RLG  INS,  GPS,  barometric  altimeter  and  radar  altimeter  avionics 
used  in  this  thesis  for  MSOFE  simulations. 
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III.  Design  Methodology  and  Error  Models 


This  chapter  describes  the  set-up  of  the  MSOFE  computer  simulation  for  the 
Landing  System  Model  (LSM)  error  model.  This  chapter  also  describes  the  technique 
used  to  determine  which  “real-world”  (true)  satellite  vehicle  (SV)  ephemeris  data  was 
used  during  MSOFE  simulation.  A  brief  description  of  the  use  of  PROFGEN  [47]  to 
generate  a  transport  flight  profile  will  also  be  discussed. 

3.1  Introduction  to  MSOFE 

The  name  "MSOFE"  is  an  acronym  meaning  "Multimode  Simulation  for  Optimal 
Filter  Evaluation."  MSOFE  is  a  general-purpose,  multimode  simulation  program  for 
designing  integrated  systems  that  employ  optimal  (Kalman)  filtering  techniques  and  for 
evaluating  their  performance  [48],  Its  general-purpose  construction  allows  specific  user 
problems  to  be  simulated  more  quickly  and  at  less  cost  than  without  its  use.  MSOFE  has 
been  designed  to  support  a  wide  variety  of  system  simulation  and  filter  evaluation  efforts. 
It  provides  two  major  operating  modes: 

1)  Monte  Carlo  simulation:  to  generate  multiple  sample  time  histories  of 
system  truth  states,  filter  states,  and  filter  estimation  errors,  including  nonlinear 
effects;  usable  for  linear  and  extended  Kalman  filters; 

2)  Covariance  simulation:  to  generate  time  histories  of  the  second-order 
statistics  (covariances)  of  system  truth  states,  filter  states,  and  filter  estimation 
errors,  under  the  assumption  of  linear  (or  linearized)  models. 

The  Monte  Carlo  and  covariance  simulation  modes  of  MSOFE  are  complementary  to  one 
another.  The  covariance  mode  can  generate  filter  performance  statistics  via  a  single  run, 
whereas  the  Monte  Carlo  mode  requires  several  sample  runs  (say,  15  or  more)  to  generate 
meaningful  statistics  for  a  given  scenario.  However,  the  covariance  mode  is  limited  to 
linear  (or  linearized)  systems,  whereas  the  Monte  Carlo  mode  can  represent  nonlinear  as 
well  as  linear  dynamic  and  measurement  processes.  In  addition,  the  Monte  Carlo  mode 
provides  better  visibility  into  the  detailed  workings  of  the  filter  models  and  computation 
processes,  and  can  easily  be  reduced  to  a  deterministic  mode  (by  suppressing  noise 
sampling)  when  required.  Monte  Carlo  runs  (15  runs)  were  solely  performed  in  this 
thesis  for  each  case  unless  otherwise  noted. 
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MSOFE  provides  a  general-purpose  simulation  environment  in  which  the  user 
embeds  a  specific  problem  by  supplying  up  to  14  problem-specific  subroutines.  The 
collective  set  of  14  user  routines  is  named  USOFE.  The  name  MSOFE  generally 
references  the  whole  program,  that  is  the  union  of  63  nonvarying  routines  in  the  CSOFE 
"core  part"  with  the  14  routines  in  the  USOFE  "user  part".  From  one  problem  to  the  next, 
the  14  routines  of  USOFE  vary  greatly,  whereas  the  63  routines  of  CSOFE  vary  only  in 
the  sizes  assigned  to  the  vectors  and  arrays.  With  MSOFE,  users  can  quickly  apply  their 
engineering  skills  to  important  filter  design  issues,  rather  than  to  the  time-consuming 
development  of  support  tools  [48]. 

The  multimode  simulation  program  MSOFE  is  part  of  an  existing  set  of  tools 
developed  by  [48]  to  support  the  design,  analysis,  and  evaluation  of  a  wide  variety  of 
integrated  systems.  Other  tools  that  were  used  in  this  thesis  were: 

•  PROFGEN  —  a  trajectory  generator  for  simulating  the  translational  and  rotational 
dynamics  of  an  aircraft  in  flight  [47], 

•  MPLOT  —  a  postprocessor  program  for  satisfying  the  plotting  needs  of  both  MSOFE 
and  PROFGEN,  e.g.  for  computing  ensemble  statistics  from  Monte  Carlo  runs  [63], 

3.1.1  MSOFE  Computer  Requirements 

MSOFE  is  written  in  ANSI  standard  FORTRAN  77  to  provide  full  portability 
across  a  wide  variety  of  computers  and  compilers.  MSOFE  was  run  on  a  486  PC  and  also 
on  a  SPARC  10  UNIX  machine  for  this  thesis.  MSOFE  is  fully  compliant  with  the  ANSI 
standard  except  for  the  way  it  manages  global  common  blocks.  These  blocks,  which  are 
constructed  in  the  form  of  individual  files,  one  file  per  block,  are  inserted  in  the  code  at 
designated  locations.  This  is  called  an  "INCLUDE"  approach,  borrowing  this  name  from 
FORTRAN  90,  where  this  ability  is  an  integral  part  of  that  standard.  The  INCLUDE 
ability  automates  common  block  array  sizing  and  aids  program  maintenance.  It  can  be 
easily  eliminated  or  modified  for  non-supporting  computers.  This  exception  to  ANSI 
compliance  was  permitted  because  of  its  usefulness  and  its  wide  availability  as  an 
extension  to  most  FORTRAN  implementations. 

The  principal  system  requirements  necessary  to  run  MSOFE  on  a  given  computer 
are: 

•  FORTRAN-77  compiler  and  linker; 

•  Ten  input/output  files  open  concurrently; 
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•  Program  and  data  memory  to  load  the  entire  program  at  one  time  (it  is  not 
overlaid):  approximately  200,000  words; 

•  Output  data  storage  of  approximately  40  (Mb)  per  hour  in  Monte  Carlo  mode 
was  not  uncommon  during  simulations  of  this  thesis  (15  Monte  Carlo  runs). 

MSOFE  allows  models  of  any  size,  limited  only  by  the  amount  of  computer 
memory  available  and  perhaps  by  array-size  restrictions  present  in  the  FORTRAN 
compiler.  There  are  no  size  restrictions  whatsoever  within  the  core  code. 

In  order  for  the  reader  to  see  the  "Big-MSOFE-Picture",  Figure  3-1  and  Figure  3-2 
(later  in  Section  3.4)  illustrate  the  overall  goal:  GPS  and  radar  altimeter  measurement 
information  must  be  fed  into  an  extended  Kalman  filter  to  determine  the  errors,  hx ,  in  the 
INS.  As  stated  earlier,  our  extended  Kalman  filter  estimates  the  true  error,  Sx ,  of  the  INS 
with  an  output  we  note  as  “Sic  ”.  Once  the  best  estimates,  8x ,  are  determined  by  our 
extended  Kalman  filter,  we  then  subtract  them  (in  a  feed-forward  approach)  from  the 
output  of  the  simulated  INS  blackbox.  The  feedforward  approach  is  utilized  in  this  thesis 
due  to  current  Federal  Aviation  Administration  (FAA)  restrictions  on  providing 
"feedback"  to  the  INS.  Without  going  into  a  lengthy  technical  discussion  about  the 
differences  between  feedforward  or  feedback,  feedforward  was  chosen  because  it  is  the 
most  conservative  choice,  especially  if  one  does  not  have  faith  that  the  extended  Kalman 
filter  feedback  corrections  will  always  be  reliable.  The  authors  belief  is  that,  at  the 
present  time,  the  FAA  does  not  want  to  "lose"  the  "pure"  INS  output  during  precision 
approaches.  Because  of  this  FAA  requirement,  a  feedforward  approach  will  be  utilized  in 
this  thesis. 


3.2  Introduction  to  PROFGEN 

PROFGEN  computes  position,  velocity,  acceleration,  attitude  and  attitude  rate  for 
an  aircraft  moving  over  the  earth  [47].  Position  is  given  as  (geographic)  latitude, 
longitude  and  altitude.  Velocity  with  respect  to  earth  is  coordinatized  and  presented  in  a 
local- vertical  frame.  Acceleration  consists  of  velocity  rates-of-change  summed  with 
Coriolis  effects  and  gravity.  Attitude  consists  of  roll,  pitch  and  yaw:  the  Euler  angles. 
[47], 
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Feedforward  corrections 


Figure  3  -  1 .  Overall  Landing  System  Model  (LSM)  Description 


PROFGEN  models  a  point  mass  responding  to  maneuver  commands  specified  by 
the  user.  These  maneuvers  are  available: 

•  vertical  turns  (pitch  up  or  down) 

•  horizontal  turns  (yaw  left  or  right) 

•  sinusoidal  "jink"  heading  changes  (oscillates  left  and  right) 

•  straight  flights  (great  circle  or  rhumb  line  path) 

PROFGEN  is  used  to  create  an  extended  flight  profile  by  concatenating  a  sequence  of 
maneuvers  chosen  from  the  basic  four.  The  user  specifies  how  long  each  maneuver  shall 
last  and  thereby  divides  the  flight  profile  into  flight  segments  [47].  Up  to  fifty  flight 
segments,  may  be  strung  together  to  produce  a  varied  total  profile.  The  final  values  of  the 
variables  in  each  segment  are  passed  along  as  the  initial  values  for  the  start  of  the  next 
segment  thereby  creating  uninterrupted  time  histories  for  all  output  variables.  The 
segments  for  the  “Tanker”  flight  profile  used  in  this  thesis  are  shown  in  Table  3-1.  The 
earth  is  modeled  as  a  perfect  ellipsoid  having  values  for  eccentricity,  semimajor  axis 
length,  spin  velocity  and  gravitational  constant  equal  to  those  of  the  DOD  World 
Geodetic  System  1972  [47]. 

A  PROF_IN  file  (PROF_IN  is  the  input  file  for  PROFGEN  [47])  was  created  for  a 
Boeing  707  (KC-135)  aircraft  based  on  inputs  from  [4,59,46].  The  actual  PROF_IN  file 
used  with  PROFGEN  can  be  found  in  Appendix  C. 
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PROFGEN  reads  in  PROF_IN,  and  outputs  a  binary  flight  file  called  "FLIGHT". 

FLIGHT  contains  the  flight  profile  variables  shown  in  Table  3-2.  Looking  at  Table  3-2 
may  seem  cryptic,  but  note  in  the  far  right  column  labeled  "WRITTEN  TO  FLIGHT"  are 
the  parameters  written  to  the  Tanker  profile  ("FLIGHT")  which  are  followed  by  the  word 
"YES";  otherwise  the  parameters  are  neither  computed  nor  written  to  the  Tanker 
"FLIGHT"  file. 

3.3  The  LSM  Computer  Model 

The  LSM  computer  model  was  derived  directly  from  the  NRS  model  [45].  Only 
the  LSM  will  be  discussed  in  detail  in  this  thesis.  For  more  information  on  the  NRS 
model,  see  [45],  The  computer  modeling  of  the  LSM  is  divided  into  two  portions,  the 
truth  model  and  the  filter  model.  The  truth  model  represents  computer-generated 
simulation  of  error  characteristics  of  avionics  "black-boxes"  and  environment,  that  would 
normally  be  found  in  the  real  world.  Because  the  information  and  data  presented  in  this 
thesis  was  accomplished  entirely  through  computer  simulation,  the  truth  model  will 
simulate  the  errors  in  true  avionics  hardware  (INS,  GPS,  Baro,  Radar  Altimeter)  black¬ 
boxes.  The  truth  model  generates  the  measurement  updates  for  the  LSM  filter,  the  true 
flight  profile  of  the  aircraft  and  a  state  variable  baseline  for  evaluating  filter  performance 
[45].  The  truth  model  consists  of  69  error  states  about  their  nominal  values.  The  filter 
model  represents  the  LSM  in  its  functional  form,  which  is  the  basis  of  the  filter  which 
could  be  hosted  on-board  an  aircraft  computer.  The  LSM  filter  model  is  a  13-state 
extended  Kalman  filter  that  has  been  developed  through  order  reduction  of  the  95-state 
truth  model  of  [45].  The  author's  approach  was  to  begin  the  filter  state  building  using  the 
fewest  possible  states  that  would  meet  precision  approach  landing  requirements.  The  13- 
state  LSM  was  chosen  as  a  first-cut  model.  An  advantage  of  using  only  13  states  is  that  it 
is  not  over-burdening  to  current  state-of-the-art  aircraft  host  computers,  keeping 
practicality  and  dollar  affordability  in  mind.  At  the  completion  of  this  thesis,  a 
determination  will  be  made  as  to  the  "usefulness"  of  the  13-state  filter.  The  tradeoffs 
involved  in  using  additional  error-states  to  increase  performance  will  be  discussed  in 
Chapter  5. 


3-5 


Segment 

Number 

Start 

Time 

(sec) 

Segment 

Length 

(sec) 

Nominal 

Path 

Centrifugal 

Accel 

(max  in  g's) 

Accel 

along 

Velocity 

vector 

ARoll 

(deg) 

APitch 

(deg) 

AHeading 

(deg) 

1 

0 

3 

STRT 

0 

0 

0 

0 

0 

2 

3 

30 

STRT 

0 

0.262 

0 

0 

0 

3 

33 

160 

VERT 

0.21 

3.80E-02 

0 

5 

0 

4 

193 

329 

STRT 

0 

1.00E-02 

0 

0 

0 

5 

522 

5 

VERT 

0.9 

-8.75E-02 

0 

-5 

0 

6 

527 

25 

HORZ 

0.9 

1  00E-02 

0 

0 

45 

7 

552 

371 

STRT 

0 

1.00E-02 

0 

0 

0 

8 

923 

35 

HORZ 

0.9 

0 

1 

0 

45 

9 

958 

400 

STRT 

0 

8.00E-03 

0 

0 

0 

10 

1358 

35 

HORZ 

0.9 

0 

1 

0 

45 

11 

1393 

340 

STRT 

0 

0 

0 

0 

0 

12 

1733 

35 

HORZ 

0.9 

0 

0 

0 

45 

13 

1768 

600 

JINK 

3.41E-02 

0 

0 

200  (sec) 

for  JINK 

2 

14 

2368 

35 

HORZ 

0.9 

0 

0 

0 

45 

15 

2403 

120 

VERT 

0.311 

0 

0 

-5 

0 

16 

2523 

35 

VERT 

0.9 

0 

0 

5 

0 

17 

2558 

372 

STRT 

0 

3.00E-03 

0 

0 

0 

18 

2930 

70 

HORZ 

0.9 

-0.125 

0 

0 

19 

3000 

40 

VERT 

0.63 

-8.50E-02 

0 

-5 

0 

20 

3040 

35 

VERT 

0.9 

0 

0 

5 

0 

21 

3075 

360 

VERT 

0.9 

1.50E-02 

0 

-1.75 

0 

22 

3435 

45 

HORZ 

0.9 

-5.00E-02 

1 

0 

-45 

23 

3480 

104 

VERT 

0.9 

0 

0 

1.75 

0 

24 

3584 

35 

HORZ 

0.9 

-5.00E-02 

1 

0 

45 

25 

3619 

39 

VERT 

0.24 

-0.1882 

0 

-3 

0 

26 

3658 

25 

HORZ 

0.9 

0 

0 

0 

45 

27 

3683 

242 

STRT 

0 

0 

0 

0 

0 

Table  3-1.  PROFGEN  Segments  for  "Tanker"  Flight  Profile 


3-6 


VARIABLE 

DIM 

PRINTED 

PRINTED 

WRITTEN 

NAME ( s ) 

TO  PROF_OUT 

TO  FLIGHT 

0 

time 

i 

TIME 

YES 

YES 

1 

terrestrial  longitude 

i 

TLON 

YES 

YES 

2 

geographic  latitude 

i 

GLAT 

YES 

YES 

3 

altitude 

i 

ALT 

YES 

YES 

4 

celestial  longitude 

i 

CLON 

5 

wander  angle 

i 

ALPHA 

YES 

YES 

6 

heading 

i 

HEAD 

7 

roll 

i 

ROLL 

YES 

YES 

8 

pitch 

i 

PITCH 

YES 

YES 

9 

yaw 

i 

YAW 

YES 

YES 

10 

terrestrial  longitude  rate 

i 

DTLON 

11 

geographic  latitude  rate 

i 

DGLAT 

12 

altitude  rate 

i 

DALT 

13 

celestial  longitude  rate 

i 

DCLON 

14 

wander  angle  rate 

i 

DALPHA 

YES 

YES 

15 

heading  rate 

i 

DHEAD 

16 

roll  rate 

i 

DROLL 

YES 

YES 

17 

pitch  rate 

i 

DPITCH 

YES 

YES 

18 

yaw  rate 

i 

DYAW 

YES 

YES 

19 

signed  earth  velocity  magnitude 

i 

VET 

20 

signed  earth  velocity  magnitude  rate 

i 

DVET 

21 

position  in  frame  i  (2) 

(3! 

3 

RI 

22 

earth  velocity  in  frame  i 

3 

VEI 

23 

inertial  velocity  in  frame  i 

(4) 

3 

VII 

24 

gravitation  in  frame  i 

3 

GNI 

25 

specific  force  in  frame  i 

(5) 

3 

FII 

26 

angular  rate,  b/i  in  frame  i 

3 

WBII 

27 

DCM  to  inertial  from  body 

3X3 

CIB 

28 

DCM  to  inertial  from  earth 

(6) 

3X3 

CIE 

29 

angular  rate,  e/i  in  frame  e 

3 

WEIE 

30 

position  in  frame  e 

3 

RE 

31 

earth  velocity  in  frame  e 

3 

VEE 

32 

inertial  velocity  in  frame  e 

3 

VIE 

33 

gravity  in  frame  e 

3 

GYE 

34 

specific  force  in  frame  e 

3 

FIE 

35 

angular  rate,  b/i  in  frame  e 

3 

WBIE 

36 

DCM  to  earth  from  body 

3X3 

CEB 

37 

DCM  to  earth  from  nav 

3X3 

CEN 

38 

angular  rate,  n/e  in  frame  n 

3 

WNEN 

39 

position  in  frame  n 

3 

RN 

40 

earth  velocity  in  frame  n 

3 

VEN 

YES 

YES 

41 

inertial  velocity  in  frame  n 

3 

VIN 

42 

gravity  in  frame  n 

3 

GYN 

43 

specific  force  in  frame  n 

3 

FIN 

YES 

YES 

44 

angular  rate,  b/i  in  frame  n 

3 

WBIN 

45 

DCM  to  nav  from  body 

3X3 

CNB 

46 

angular  rate,  b/n  in  frame  b 

3 

WBNB 

47 

position  in  frame  b 

3 

RB 

48 

earth  velocity  in  frame  b 

3 

VEB 

49 

inertial  velocity  in  frame  b 

3 

VIB 

50 

gravity  in  frame  b 

3 

GYB 

51 

specific  force  in  frame  b 

3 

FIB 

52 

angular  rate,  b/i  in  frame  b 

3 

WBIB 

Table  3-2.  PROFGEN  Right  Profile  Outputs 
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X 


Figure  3-2.  Truth  and  Filter  Model  Block  Diagram 


The  block  diagram.  Figure  3-2  explains  how  the  filter  and  truth  models  interact 
in  the  MSOFE  computer  simulation.  PROFGEN  [47]  provides  a  simulated  flight  profile 
and  the  U.S.  Coast  Guard  GPS  Bulletin  Board  Service  (BBS)  [24]  provides  true  SV 
ephemeris  data  for  any  SV.  Use  of  the  "real-world"  ephemeris  replaced  the  prior 
FORTRAN  ORBIT  functions  used  by  past  researchers  at  AFIT  [45,58,3].  The  best  four 
SV  were  chosen  by  using  System  Effectiveness  Model  (SEM)  software  [18]  and  selecting 

the  best  (lowest)  position  dilution  of  precision  (PDOP).  With  this  information,  the  truth 
model  is  able  to  simulate  a  real  world  INS  navigation  solution,  x+bxINS ,  and  generate 

the  real  world  GPS  and  radar  altimeter  measurements,  RGPS  and  Rrait  respectively.  The 
LSM  filter  in  Figure  3  -  2  is  represented  by  the  Kalman  filter  block.  The  corrections  from 

the  LSM  filter  are  subtracted  from  the  INS  navigation  solution  to  generate  the  best 
possible  navigation  solution  available,  x  =  x+bxINS-dxINS  [45].  The  switch  in  Figure  3 

-  2  does  not  imply  "either/or",  instead  it  implies  use  of  radar  altimeter  measurements  as 
well  as  GPS  receiver  outputs  can  be  used.  Now  that  the  MSOFE  implementation  of  the 
LSM  filter  has  been  explained,  the  truth  and  filter  models  for  the  GPS,  radar  altimeter  and 
the  INS  subsystems  will  be  described. 
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3.4  LSM  Model  Description 


This  section  presents  the  truth  and  filter  model  propagation  and  measurement 
equations,  (2.35)  and  (2.36),  respectively.  The  following  presentation  will  be  divided  up 
by  navigation  subsystems  with  most  parts  taken  directly  from  [45].  First  the  INS  portion 
of  the  equations  will  be  presented,  then  the  GPS,  followed  by  the  radar  altimeter. 

Before  the  different  navigation  subsystems  are  individually  described,  the  high-level  state 
and  measurement  equations  are  provided  for  the  LSM  filter,  followed  by  those  for  the 
truth  model.  Equations  (3.1)  and  (3.2)  show  how  the  different  navigation  subsystems 
models  combine  to  form  a  single  LSM  filter  model: 


FINSf 

0 


0 

FGPSf 


+ 


wINSf 

wGPSf 


HlNSf 

HGPSf 


&/  + 


vINSf 
v  INS j 


(3.1) 


(3.2) 


As  stated  earlier,  the  overall  filter  model  consists  of  13  states;  1 1  for  the  INS  and  2  for  the 
GPS.  Table  A.4  in  Appendix  A  provides  a  description  of  the  13-state  vector,  bxj: , 

implemented  in  the  filter  model.  References  to  further  descriptions  of  the  sub-matrices  in 
the  filter  equations  can  be  found  in  Table  3-3.  It  should  be  noted  that  the  barometric 
altimeter  aiding  measurements  are  considered  to  be  INS  measurements,  while  the  GPS 
and  radar  altimeter  measurements  are  the  respective  updates  for  the  baro/inertial  system 
from  the  GPS  and  radar  altimeter. 


The  propagation  and  measurement  equations  for  the  LSM  truth  model  is 
presented  in  similar  fashion  below: 


1 _ 

FINStl 

0 

W  Filter 

8? 

II 

0 

finsi2 

0 

bxt  + 

WINS, 

0 

0 

fgps,  _ 

_WGPS, 
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Filter  Model 

Location  of  Description 

Truth  Model 

Location  of  Description 

FINSf 

Section  3.4. 1.3 

F INS  Filter 

Section  3.4. 1.3,  3.4.2. 1 

- 

- 

FINSn 

Section  3.4. 1.2 

- 

- 

F1NS,2 

Section  3.4. 1.2 

FGPSf 

Section  3.4.2.2 

fgps, 

Section  3.4.2. 1 

wINSf 

Section  3.4. 1.3 

WINS, 

Section  3.4. 1.2 

wGPSf 

Section  3. 4. 2.2 

WGPS, 

Section  3.4.2. 1 

HlNSf 

Section  3.4. 1.4 

Hjns, 

Section  3.4. 1.4 

HGPSf 

Section  3.4. 2.3 

hgps, 

Section  3. 4.2.3 

Table  3-3.  References  for  the  Sub-Matrices  of  the  LSM  Truth  and  Filter 


6zf  = 


hins, 

8x,  + 

VINS, 

hgps, 

l 

VGPS, 

(3.4) 


The  LSM  truth  model  consists  of  the  original  thirteen  states  of  the  filter  model 
(represented  by  FFl!ter  and  wmer ),  augmented  by  additional  INS  and  GPS  states  (the 
radar  altimeter  measurements  were  modeled  as  corrupted  only  by  white  noise,  thus  no 
additional  states  were  necessary).  The  total  number  of  states  for  the  navigation 
subsystems  is  69;  39  INS  states,  and  30  GPS  states.  Tables  A.1-A.3,  in  Appendix  A  gives 
a  full  description  of  each  individual  state  of  the  truth  model.  Also  Tables  B.1-B.4  and 
Tables  B.6-B.7  in  Appendix  B  have  a  complete  listing  of  the  components  of  the  F  and  the 
Q  noise  strengths  associated  with  the  w  vector  components  in  Equation  (3.3). 

While  the  first  thirteen  states  of  the  filter  model  are  nearly  identical  to  the  first 
thirteen  states  of  the  truth  model,  there  is  one  crucial  difference  [45].  The  filter  model 
dynamics  driving  noise  and  measurement  noise  do  not  correlate  exactly  with  those  of  the 
first  thirteen  states  of  the  truth  model.  The  filter  model  noise  statistics  values  have  been 
altered  to  achieve  good  tuning  against  the  truth  model  [40],  The  following  sections  will 
provide  a  detailed  presentation  into  the  exact  make-up  of  the  truth  and  filter  model 
propagation  and  measurement  equations  for  all  navigation  subsystems  used  in  this  thesis. 
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3.4.1  The  Inertial  Navigation  System  (INS)  Model. 


This  section  presents  the  truth  and  filter  models  used  for  the  INS.  The  INS  model 
is  a  strapped-down  wander  azimuth  system  that  senses  aircraft  motion  via  gyros  and 
accelerometers  and  is  used  as  the  primary  source  for  navigation  [45].  The  INS  model  has 
been  derived  from  a  medium  accuracy  RLG  INS  93-state  model  [1,38].  First,  the  original 
93-state  model  will  be  presented,  followed  by  the  reduced-ordered  39-state  truth  and  1 1- 
state  filter  models.  After  the  INS  truth  and  filter  state  equations  have  been  defined, 
barometric  altimeter  measurement  equations  will  be  presented. 

3. 4.1.1  The  93-State  LN-93  Error  Model. 


The  93-state  Litton  INS  MSOFE  computer  model  has  been  generated  by  the 
Wright  Laboratories,  Avionics  Directorate,  Avionics  System  Integration  and  Research 
Team  (ASIRT).  Their  development  uses  both  past  AFIT  research  and  INS  vendor  [38] 
documentation  to  “fine-tune”  past  modeling  efforts  [50,58,65].  The  93-state  model 
generates  a  high  number  of  documented  error  sources  that  are  found  in  the  Litton  wander- 
azimuth  LN-93  INS  [38],  These  errors  are  described  using  six  categories  of  states  [45]: 


8x  =  [Sxf  SjcJ  8x3r  8xJ  8xJ  8xJ  j7 


(3.5) 


where  8x  is  a  93  x  1  column  vector  and: 

Sxj :  represents  the  "general"  error  vector  containing  13  position,  velocity,  attitude,  and 
vertical  channel  errors  (representative  of  a  Pinson  model  of  INS  error 
characteristics). 

&c2 :  consists  of  16  gyro,  accelerometer,  and  baro-altimeter  exponentially  time- 
correlated  errors,  and  "trend"  states.  These  states  are  modeled  as  first  order 
Markov  processes  in  the  truth  (system)  model. 

8x3  :  represents  gyro  bias  errors.  These  18  states  are  modeled  as  random  constants  in 

the  truth  model. 

8x4 :  is  composed  of  the  accelerometer  bias  error  states.  These  22  states  are  modeled  in 

exactly  the  same  manner  as  the  gyro  bias  states. 

8x5 :  depicts  accelerometer  and  gyro  initial  thermal  transients.  The  6  thermal  transient 

states  are  first  order  Markov  processes  in  the  system  model. 
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S*6  :  models  the  gyro  compliance  errors.  These  18  error  states  are  modeled  as  biases  in 
the  system  model. 

The  93-State  Litton  model  state  space  differential  equation  is  given  by: 
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A  full  description  of  the  sub-matrices  for  this  equation  is  given  in  the  Litton  LN-93 
manual  [38],  This  large  state  model  represents  the  most  accurate  model  available  for  the 
LN-93  navigation  errors  [45]. 

3.4. 1.2  The  39-State  INS  Truth  Model. 


While  the  93-state  model  is  a  very  accurate  representation  of  the  INS  error 
characteristics,  the  high  dimensionality  of  the  state  equations  makes  the  model  very  CPU¬ 
intensive  for  “first-look”  projects.  The  intent  of  this  thesis  is  not  to  generate  a  high 
fidelity  model  of  the  Litton  LN-93,  but  rather  to  evaluate  performance  characteristics 
associated  with  a  particular  class  of  INS  (medium  precision  or  lower  precision).  Previous 
AFIT  theses  have  demonstrated  that  reduced-ordered  truth  models  can  be  used  in  place  of 
the  93-state  truth  model  without  losing  a  significant  degree  of  accuracy  [45,52,50]. 
Therefore  the  INS  truth  model  has  been  reduced  to  a  39-state  model.  The  reduced- 
ordered  model  retains  only  the  truly  essential  states  from  Equation  (3.6).  The  truth  model 
state  space  equation  is  defined  in  Equation  (3.7): 
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(3.7) 


It  should  be  noted  that  the  INS  truth  state  vector  8*,  is  a  39-state  vector.  The  four 
components  of  8*  do  not  directly  correlate  to  the  first  four  components  of  the  93-state 
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Litton  model  [45].  For  a  complete  listing  of  the  39  states  and  how  they  relate  to  those  in 
[38],  see  Tables  A.  1  and  A.2  in  Appendix  A. 

3.4. 1.3  The  11-State  INS  Filter  Model. 

The  INS  filter  model  retains  the  essential  states  from  the  39-state  truth  model. 
Through  past  AFIT  research,  the  1 1 -state  INS  filter  has  been  shown  to  perform 
adequately  when  given  frequent  GPS  measurement  updates  [45,50,65].  Table  A.4  in 
Appendix  A  shows  the  1 1  states  used  for  the  INS  filter  model.  The  final  INS  filter 
dynamics  submatrix,  F,  as  well  as  process  noise  strength  Q  and  measurement  noise 
covariance  R,  can  be  found  in  Appendix  B. 

3.4. 1.4  INS  Measurement  Model. 

The  two  measurements  that  are  used  to  update  the  filter  are  the  barometric 
altimeter  and  the  radar  altimeter.  The  barometric  altimeter  signal  is  used  to  correct  for 
inherent  instabilities  of  the  vertical  channel  in  the  filter.  The  radar  altimeter  is  used 
during  landing  approaches  when  altitudes  are  below  3000  feet  above  ground  level  (AGL). 
First  the  barometric  altimeter  measurement  will  be  presented,  followed  by  the  radar 
altimeter  measurement.  It  should  be  noted  that,  since  the  LSM  filter  is  an  error  state  filter, 
it  is  necessary  to  develop  difference  measurement  update  equations  for  all  the 
measurements.  The  barometric  altimeter  measurement  equation  is  based  on  the 
difference  between  the  INS-predicted  altitude,  Altm  and  the  barometric  altimeter- 
predicted  altitude  AltBar: 

5 ZAlt  =  Alt  INS  -  AltBar  @.8) 

Therefore  it  is  necessary  to  develop  the  two  separate  measurement  signals  that  will  be 
differenced  to  attain  the  proper  measurement  update  for  the  error  state  filter  [45].  The 
INS-predicted  altitude  is  the  sum  of  the  true  altitude,  ht,  and  the  INS  error  in  vehicle 
altitude  above  the  reference  ellipsoid,  5 h .  The  barometric  altimeter  reading  is  modeled 
as  the  sum  of  the  true  altitude,  ht,  the  total  time-correlated  error  in  the  barometric 
altimeter,  dhB ,  and  a  random  measurement  noise,  v.  The  difference  measurement  update 

signal  is  formed  in  Equation  (3.9)  by  subtracting  the  INS-predicted  altitude  from  the 
barometric  altimeter  altitude: 


3-13 


^zAlt  -  Ah INS  ~  Alt  Bar 


=  [ht  +6h\-[ht+  bhB  -  v]  (3.9) 

=  5  h-bhB  +  v 

Note  that,  since  v  is  assumed  zero  mean,  white  gaussian  noise,  statistically  speaking,  one 
can  choose  v  with  either  a  plus  "+"  or  minus  coefficent .  The  author  chooses  the 
coefficent  carefully  so  that  the  end  result  shows  a  plus  "+v"  sign. 

This  completes  the  presentation  of  the  INS  truth  and  filter  state  equations  as  well 
as  the  INS  measurement  equations.  The  next  section  will  develop  similar  equations  for 
the  radar  altimeter  used  in  this  thesis. 

3.4.2  Radar  Altimeter  Model. 

As  a  “first-cut”  model  of  the  radar  altimeter,  the  measurement  equation  is  based 
on  the  difference  between  the  INS  predicted  altitude,  AIt1NS  and  the  radar  altimeter 
predicted  altitude,  AltRait: 

§ZALT  ~  ALT,ns  -  ALTRalt 

=  [h,+8h]-[ht-v]  (3.10) 

-bh  +  v 

Note  that  the  errors  in  the  radar  altimeter  are  represented  totally  as  white  noise,  with  no 
time-correlated  component  at  all.  Though  admittedly  only  a  first-cut  model,  it  should  be 
sufficient  to  demonstrate  important  performance  trends. 

The  radar  altimeter  measurement  noise  variance,  RFi[ter  or  Rjrue,  is  a  function  of 
aircraft  altitude  above  ground  level  (AGL).  The  filter  model  noise  variance  from  [28]: 

RTruth  =  {[0.01]2  * [ Radar  AltitudeTrue{AGL)}2  }  +  0. 25 Bias  (3.11) 

and  the  truth  model  noise  variance  is  the  same: 

RTruth  =  {[0.01]2  * [ Radar  AltitudeTrue{AGL) ]2 }  +  0. 25 Bias  (3.12) 

Note  that  Rputer  and  Rjruth  are  both  time- varying  rather  than  constant,  due  to  the  altitude 
dependency. 
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This  completes  the  presentation  of  the  radar  altimeter  measurement  and  noise 
variance  equations.  The  next  section  will  develop  truth,  filter  and  measurement  equations 
for  the  Global  Positioning  System  (GPS)  used  in  this  thesis. 

3.4.3  The  Global  Positioning  System  ( GPS)  Model. 

The  GPS  navigation  system  used  is  based  on  electromagnetic  signals  transmitted 
from  orbiting  GPS  satellites.  This  model  has  been  developed  throughout  research  at 
AWT,  and  many  of  its  fundamental  concepts  are  addressed  in  a  variety  of  sources 
[39,50,56,58].  GPS  generates  navigation  information  by  acquiring  the  range  to  multiple 
satellites  of  known  position,  called  “ pseudoranges ”.  The  navigation  information  passed 
to  the  LSM  filter  is  the  respective  range  and  ephemeris  data  position  to  each  of  four 
satellites  [45,39].  The  next  three  sections  will  present  all  the  necessary  equations  to 
define  the  GPS  truth  and  filter  error  models  fully. 

3.4.3. 1  The  30-State  GPS  Truth  Model 

There  are  five  types  of  error  sources  that  are  modeled  in  the  GPS  truth  model  state 
equations.  The  first  two  states  represent  the  errors  in  the  user  clock  and  are  modeled  as 
follows: 


where 
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range  equivalent  of  user  clock  bias 
velocity  equivalent  of  user  clock  drift 


(3.13) 


The  initial  state  estimates  and  covariances  for  these  states  were  chosen  to  be  consistent 
with  previous  AWT  research,  [45,50,56,58]  and  are: 
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Because  these  error  sources  are  a  function  of  the  user  equipment,  they  are  common  to  all 
the  satellite  vehicles.  The  remaining  five  sources  of  errors  are  unique  to  each  satellite 

vehicle  (SV),  based  on  their  individual  equipment  and  their  position  with  respect  to  the 
user.  The  first  SV-specific  error  source  is  the  code  loop  error,  8/?c/oop .  Although  the 

code  loop  is  part  of  the  user  equipment  shared  by  all  the  SV's,  its  error  magnitude  is 

relative  to  each  SV.  The  second  and  third  SV-specific  errors  are  the  atmospheric 
interference  with  the  EM  signals,  &Rion  and  hRtrop ,  as  related  to  the  ionospheric  and 

tropospheric  delay  in  the  signal's  propagation.  The  code  loop  error,  tropospheric  delay, 
and  ionospheric  delay  are  all  modeled  as  first  order  Markov  processes  with  time  constants 
shown  in  Equation  (3.16),  consistent  with  previous  AFIT  research  [58,65,50,45].  All 
three  are  driven  by  zero-mean  white  Gaussian  noise  with  strengths  shown  in  Equation 

(3.19).  The  fourth  SV-specific  error  source  is  due  to  inaccuracies  in  the  clocks  on  board 
the  SV's,  5 Rsdk  and  the  final  error  source  is  based  on  line-of-sight  errors  between  the 

SV's  and  the  receiver,  8x5  ,  8ys  ,  hzs. ,  respectively.  The  model  for  these  states  shown  in 
Equation  (3.16): 

8Rcioop)  r  -1  0  0  0  0 

SRtrop  0  “500  0  0  0 

S^ion  0  0  “1500  0  0 

'  8i?ScW.  1=  0  0  0  0  0 

8i.  0  0  0  0  0 

■*» 

dys.  0  0  0  0  0 

dL  0  0  0  0  0 

**/  L 

where  the  initial  covariance  for  the  states  is  given  by: 

0.25  ft2  0  0  0  0  0 

0  1.0  ft2  0  0  0  0 

0  0  1.0  ft2  0  0  0 

Pgps=  0  0  0  0.25  ft2  0  0 

0  0  0  0  25ft2  0 
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0  (3.17) 
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and  mean  values  and  strengths  of  the  dynamics  driving  noise  are  given  by: 
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A  quick  reference  of  the  truth  model  non- zero  GPS  dynamics  matrix  components  is 
provided  in  Tables  B.4  of  Appendix  B.  This  ends  the  description  of  the  30-state  truth 
model.  Now  the  filter  model  will  be  presented. 


3A.3.2  The  2-State  GPS  Filter  Model 

Various  research  efforts  have  shown  that  two  states  provide  a  sufficient  model  for 
GPS  [50,52,45].  The  primary  argument  is  that  the  errors  modeled  by  the  other  28  states 
are  small  when  compared  to  the  two  states  common  to  all  SV's.  By  adding  dynamics 
driving  noise,  of  strength  Q,  and  re-tuning  the  filter,  the  overall  performance  of  the  LSM 
can  be  maintained  with  the  significantly  reduced-order  model  of  Equation  (3.20): 

5/W 1  = 


0  1 
0  0 


SR, 


■elk,. 


m 


clk„ 


I WK, 

\WR, 


(3.20) 


The  values  implemented  for  the  dynamics  driving  noise  strengths  can  be  found  in  Table 
B.8  of  Appendix  B.  It  should  be  noted  that,  in  the  tuning  process,  the  measurement  noise 
covariance  values  R  (as  shown  in  Appendix  B)  have  also  been  adjusted  to  achieve 
adequate  tuning  of  the  filter  [40].  This  completes  the  description  of  the  GPS  filter  model. 
The  next  section  presents  the  GPS  measurement  equations  for  both  the  truth  and  the  filter 
models. 
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3.43.3 


GPS  Measurement  Model 


There  are  four  GPS  measurement  updates,  one  from  each  of  the  satellite  range 
signals  received  by  the  LSM  filter.  These  measurement  updates  are  once  again  difference 
measurements.  First  the  GPS  truth  model  difference  measurement  will  be  fully 

presented,  followed  by  a  brief  description  of  the  filter  measurement.  The  GPS  difference 
measurement  is  formed  by  taking  the  difference  of  the  INS -calculated  pseudorange,  RlNS 

and  actual  pseudorange,  RGPS . 


&zgps  ~  rins  ~  rgps  (3-21) 

The  real  pseudorange,  RGPS  is  the  sum  of  the  true  range  form  the  user  to  the  satellite  plus 
all  the  errors  in  the  pseudorange  signal  propagation. 

RC,PS  =  Rt  +  &Rcloop  +  &Rtrop  +  &Rion  +  &RSclk  +  &RUclk  ~  v  (3  -22) 


where 


rgps  - 


Rt  = 

^ Rcloop  ~ 
Rtrop  ~ 


8 


5RSclk  ~ 

&RUctk  ~ 


V 


GPS  pseudorange  measurement,  from  SV  to  user 

true  range,  from  SV  to  user 

range  error  due  to  code  loop  error 

range  error  due  to  tropospheric  delay 

range  error  due  to  ionospheric  delay 

range  error  due  to  SV  clock  error 

range  error  due  to  use  clock  error 

zero-mean  white  Gaussian  measurement  noise 


Note  that,  since  v  is  assumed  zero  mean,  white  gaussian  noise,  statistically  speaking,  one 
can  choose  v  with  either  a  plus  "+"  or  minus  coefficent .  As  in  equation  (3.9),  the 
author  chooses  the  coefficent  carefully  so  that  the  end  result  shows  a  plus  "+v"  sign  in 
equation  (3.25)  and  equation  (3.26). 

The  second  source  of  a  range  measurement  is  the  INS  itself,  R/NS  [45].  R;ns  is  the 
difference  between  the  LSM-calculated  position,  Xv  and  the  satellite  position  from  the 
ephemeris  data  Xs.  This  difference  vector  is  represented  below  in  the  ECEF  frame: 
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(3.23) 
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An  equivalent  form  for  Equation  (3.23)  is 

rins  =  ^j(xu  ~xs)2  +(yu  ~ys)2  +(zu  ~zs)2  (3.24) 


Based  on  Assumption  12  from  Chapter  1,  Equation  (3.24)  can  be  approximated  and 
rewritten  in  terms  of  the  true  range  and  a  truncated  first-order  Taylor  series,  with 
perturbations  representing  the  errors  in  Xv  and  Xs: 


rins  -Rt  + 


dR!NS(Xs,  X;;) 


dX, 


u 


(Xs.Xu\ 


■bxu  + 


X(j) 


dXc 


(XsXu)ri 


■SXs  (3.25) 


The  solution  for  R,NS  is  found  by  substituting  Equations  (3.24)  into  Equation  (3.25)  and 
evaluating  the  partial  derivatives  to  get  [45]: 
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(3.26) 


Finally,  the  GPS  pseudorange  truth  model  difference  measurement  is  given  as: 
&zgps,  =  rins  ~  rgps 
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(3.27) 
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The  user  position  errors  in  Equation  (3.27)  can  be  derived  from  the  first  three  states  of  the 
filter  or  truth  model  using  an  orthogonal  transformation  [9,45]. 

The  filter  model  for  the  GPS  measurement  will  now  be  derived.  Since  the  filter 
model  does  not  contain  the  states  for  the  errors  in  the  satellite  position,  these  terms  are 
removed  from  the  equation.  The  filter  model  measurement  equation  can  therefore  be 
written  as: 


$ ZGPSf  ~  R/NS  ~  RGPS 


xs~xu 
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(3.28) 


The  filter  measurement  noise  variance,  R,  will  be  tuned  to  attain  adequate  performance 
despite  the  reduction  in  order  from  the  truth  model  and  the  Taylor  series  approximation. 
The  measurement  noise  variances  for  both  the  filter  and  the  truth  model  equations  are 
provided  in  Table  B.9  of  Appendix  B.  This  completes  the  description  of  the  GPS 
measurement  equations  and  the  entire  LSM  filter  and  truth  model  equations. 

3.5  Chapter  Summary 

This  chapter  presents  the  modeling  set-up  of  the  LSM  MSOFE  computer 
simulation.  An  introduction  to  MSOFE  and  PROFGEN  is  provided.  The  truth  model 
and  filter  model  propagation  and  measurement  equations  are  described  for  the  INS/Baro, 
Radar  Altimeter,  and  GPS  subsystems.  A  "first-cut"  model  of  the  radar  altimeter  model 
has  been  presented.  The  radar  altimeter  measurements  should  play  a  key  role  in  allowing 
our  aircraft  to  meet  precision  approach  requirements.  The  INS/Baro  and  GPS  truth  model 
is  located  in  tabular  form  in  Appendix  A.  The  dynamic  submatrices  FFilter  FINStl,  FINSt2, 
and  Fgps,  and  process  noise  strength  and  measurement  noise  covariance  matrices  for 
filter  and  truth  models  are  presented  in  Appendix  B.  Results  and  analysis  of  the  LSM 
simulation  are  presented  in  Chapter  4. 
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IV.  Results  and  Analysis 


This  chapter  presents  results  and  analysis  of  the  following  items: 

4.1  "Tanker"  trajectory  (detailed  KC-135  flight  profile)  created  using 
PROFGEN  [47], 

4.2  SV  Ephemeris  Data,  (GDOP,  PDOP,  HDOP,  VDOP  SV  numbers)  using 
SEM  3.6  [18]. 

4.3  Development  of  The  Three  INS's 

4.4  Filter  Tuning  Example  (using  Case  I  error  states, ).  See  Table  4-1. 

4.5  Case  I- VI  performance  results  (includes  radar  altimeter  aiding). 

4.6. 1  Case  VII-IX  performance  results  (includes  single  pseudolite  aiding). 

4.6.2  Case  X-XII  performance  results  (includes  radar  altimeter  and  pseudolite 
aiding). 

4.7  Case  XIII-XV  performance  results  (includes  one  GPS  time-out  for  the 
three  INS  accuracies:  0.4  nm/hr,  2.0  nm/hr  and  4.0  nm/hr). 

4.8  CASE  XVI  performance  results  (includes  two  GPS  time-outs  with  radar 
altimeter  aiding,  during  final  approach  using  a  low  accuracy  (4.0  nm/hr) 
INS. 


As  a  helpful  aid.  Tables  4-1  and  4-2  show  Cases  I-IX  and  Cases  X-XVI  integration 
comparisons,  respectively. 


Case  I 

Case  II 

CaseHI 

Case  IV 

Case  V 

Case  VI 

Case  Vn 

Case  Vin 

Case  IX 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

0.4  nm/hr 

CEP  INS 

■K33H 

4.0  nm/hr 

CEP  INS 

4.0  nm/hr 

CEP  INS 

0.4  nm/hr 

CEP  INS 

2.0  nm/hr 

CEP  INS 

4.0  nm/hr 

CEP  INS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

Radar 

Altimeter 

Radar 

Altimeter 

— 

Radar 

Altimeter 

Pseudolite 

Pseudolite 

Pseudolite 

Table  4-1.  Case  I-IX  Integration  Comparisons 
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Case  X 

Case  XI 

Case  XII 

Case  XIII 

Case  XIV 

Case  XV 

Case  XVI 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

Barometric 

Aldmeter 

Barometric 

Altimeter 

Barometric 

Altimeter 

0.4  nm/hr 

CEP  INS 

2.0  nm/hr 

CEP  INS 

4.0  nm/hr 

CEP  INS 

0.4  nm/hr 

CEP  INS 

2.0  nm/hr 

CEP  INS 

4.0  nm/hr 

CEP  INS 

4.0  nm/hr 

CEP  INS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

P-Code 

GPS 

Radar 

Aldmeter 

and 

Pseudolite 

Radar 

Altimeter 

and 

Pseudolite 

Radar 

Aldmeter 

and 

Pseudolite 

None 

None 

None 

Radar 

Altimeter 

No  GPS 

Outage 

No  GPS 

Outage 

No  GPS 

Outage 

Single  GPS 
Outage 

Single  GPS 
Outage 

Single  GPS 
Outage 

Double 

GPS 

Outages 

Table  4-2.  Case  X-XI  Integration  Comparisons 


Tables  4-1  and  4-2  show  a  total  of  sixteen  different  cases  that  were  completed 
in  this  thesis.  All  cases  use  the  Tanker  (KC-135)  flight  profile  of  Section  4. 1  and  the  P- 
Code  receiver  always  uses  four  SVs  overhead.  Cases  I  -  VI  show  how  an  the  additional 
radar  altimeter  measurements  can  aid  three  different  INS’s  using  baro  altimeter  and  a  P- 
Code  GPS  receiver.  Cases  VII  -  IX  show  how  use  of  a  single  pseudolite  measurement  in 
close  proximity  to  a  runway  can  aid  three  different  INS's  using  baro  altimeter  and  P-Code 
GPS  aiding.  Cases  X  -  XII  show  the  performance  enhancements  of  using  a  single 
pseudolite  measurement  and  radar  altimeter  measurements.  The  results  of  Cases  X  -  XII 
can  be  compared  with  Cases  (I,  III,  V)  -  no  radar  altimeter  or  pseudolite,  Cases  (II,  IV, 

VI)  -  just  radar  altimeter,  and  Cases  (VII,  VIE,  IX)  -  just  pseudolite,  to  see  the  additional 
benefit  of  each  measurement  source.  Case  XEI  -  XV  shows  simply  a  single  GPS  outage 
with  the  0.4/2/4  nm/hr  INSs.  Use  of  the  radar  altimeter  occurs  during  final  approach. 
Case  XVI  shows  a  double  GPS  outage;  where  use  of  radar  altimeter  occurs  during  the 
aircraft  landing.  Only  the  4  nm/hr  INS  is  evaluated  for  this  special  case.  Also,  all  results 
will  be  compared  with  Table  4-3,  illustrating  landing  performance  la  requirements  at 
respective  decision  heights. 
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Precision  Approach  Parameters  (in  feet,  ail  1 -sigma  values) 


Category 

Decision  Height 

Azimuth 

Elevation 

1 

200  feet 

+/-28.1 

+/-  6.8 

II 

1 00  feet 

+/-  8.6 

+/-  2.8 

III 

50  feet 

+/-  6.8 

+/- 1.0 

Table  4-3  Precision  Approach  Accuracy  Requirements  at  Decision  Heights 


4.1  The  Tanker  (KC-135)  flight  profile 

A  PROFJN  file  (PROF_IN  is  the  input  file  for  PROFGEN  [47])  was  created  for  a 
Boeing  707  (KC-135)  aircraft  based  on  inputs  from  [4,59,46].  The  actual  PROF_IN  file 
used  with  PROFGEN  can  be  found  in  Appendix  C.  PROFGEN  reads  in  PROF_IN,  and 
outputs  a  binary  flight  file  called  "FLIGHT".  FLIGHT  contains  the  flight  profile 
variables  shown  in  Chapter  3,  Table  3-2. 

The  Tanker  flight  profile  is  simulated  to  start  at  WPAFB,  OH.  with  the  following 
initial  conditions: 


•  Runway  23  R  (Heading  232°) 

•  Initial  Latitude  39.055°  N 

•  Initial  Longitude  84.033°  W 

•  Initial  altitude  825  feet 

•  Aircraft  take-off  (rotation)  speed  150  knots 

•  Initial  wander  angle  (alpha)  45° 

As  described  in  Chapter  3,  PROFGEN  outputs  are  concatenated  segments.  See  Table  3-1 
in  Chapter  3.  Using  Table  C  -  1  of  Appendix  C  and  Table  3  -  1  of  Chapter  3,  the  Tanker 
flight  profile  will  now  be  defined: 

•  The  tanker  aircraft  rolls  down  the  runway  for  30  seconds;  during  this  time  its  velocity 
transitions  from  0  knots  to  150  knots  at  which  time  aircraft  rotation  occurs. 

•  The  tanker  climbs  at  a  +5°  pitch  over  the  next  494  seconds,  obtaining  an  altitude  of 
19,996  feet  and  velocity  of  326  knots,  at  an  elapsed  time  of  527  seconds. 

•  The  tanker  makes  a  +45°  (rolled  to  the  right)  turn  and  continues  straight,  with  a 
gradual  acceleration,  maintaining  an  altitude  of  19996  feet.  Final  velocity  is  393 
knots  at  an  elapsed  time  of  923  seconds. 
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•  The  tanker  makes  another  +45°  turn  and  continues  straight  with  a  gradual 
acceleration,  maintaining  an  altitude  of  19,996  feet.  Final  velocity  is  453  knots  at  an 
elapsed  time  1393  seconds. 

•  The  tanker  makes  another  +45°  turn  and  immediately  begins  a  600-second  segment  of 
mild  horizontal  jink  maneuvers  maintaining  an  altitude  of  19,996  feet  (the  jink 
maneuver  is  a  sinusoidal  heading  change  of  ±  2  degrees,  every  200-second  period). 
Final  velocity  is  454  knots,  at  an  elapsed  time  2368  seconds.  The  jinking  maneuver  is 
clearly  evident  in  the  Tanker  flight  profiles  of  Appendix  D:  Figure  D  -  2  (Wander 
Azimuth  Angle),  Figure  D  -  3  (X-Velocity),  and  Figure  D  -  4  (Roll  and  Heading). 

•  The  tanker  makes  another  +  45  degree  turn  and  immediately  begins  to  descend  and 
maintain  an  altitude  of  12255  feet.  Final  velocity  is  453  knots  at  an  elapsed  time  of 
2558  seconds. 

Note:  It  was  during  this  maneuver,  that  the  filter  (for  all  cases)  showed  a  degradation  in 
latitude  error  and  slight  degradation  in  longitude  error.  Possibly,  a  more  realistic  Tanker 
profile  would  have  decreased  velocity,  before  making  such  a  large  altitude  descent 
(though  the  tanker  does  not  exceed  a  lg  maneuver  at  any  given  time).  This  profile  is 
meant  to  be  realistic;  there  was  no  attempt  to  put  in  maneuvers  to  enhance  identifibility  of 
certain  states  (i.e.,  no  42-minute  "legs",  butterfly  patterns,  etc.). 

•  Deceleration  now  begins,  with  a  +90°  turn,  followed  by  a  -5°  pitch  (down)  to  an 
altitude  of  10,678  feet.  Final  velocity  is  244  knots  at  an  elapsed  time  of  3075 
seconds. 

•  Slight  acceleration  followed  by  gradual  deceleration  begins  with  a  -45°  turn  and  a 
negative  1.75°  pitch  (down)  to  an  altitude  of  4425  feet.  Final  velocity  is  304  knots  at 
an  elapsed  time  of  3584  seconds. 

•  Gradual  deceleration  continues  with  a  +45°  turn,  followed  by  a  -3°  glideslope  pitch 
(down).  True  altitude  mean  sea  level  "MSL"  is  3799  feet  or  2974  feet  above  ground 
level  "AGL",  velocity  is  131  knots  at  an  elapsed  time  of  3658  sec. 

•  The  tanker  makes  its  final  +45°  turn  to  align  with  the  WPAFB  runway  (23R  at  232° 
heading).  The  -3  degree  glideslope  continues,  until  the  tanker  aircraft  touches  down 
at  an  altitude  of  825  feet  (MSL)  or  0  feet  (AGL).  Aircraft  velocity  while  on  final 
approach  is  132  knots  at  all  times.  Touchdown  occurs  at  an  elapsed  time  of  3912.5 
sec. 

The  aircraft  reaches  Cat  I,  n,  HI  decision  heights  at  the  following  flight  profile  times 
indicated  in  Table  4-4. 


4-4 


Category  of 
Precision 
Approach 

Decision 
Height  in  feet 

Elapsed 
Time  in 
seconds 

Aircraft 
Altitude  (AGL) 
in  feet 

Aircraft 
Altitude  (MSL) 
in  feet 

I 

200 

3896 

200 

1025 

II 

100 

3904 

100 

925 

III 

50 

3908 

50 

875 

Table  4-4.  Category  I,  II  and  III  Decision  Height  Profile  Times 


Note:  Vertical  velocity  (Z-velocity)  =  -11.6975  ft/sec  =  6.92465  knots  -7.9687  mph 
during  the  final  approach.  As  stated  earlier  in  Chapter  3,  the  PROFGEN  output  file 
"FLIGHT"  is  used  as  the  truth  model  trajectory  in  the  MSOFE  LSM  simulation. 


4.2  Satellite  Vehicle  Data  Using  SEM3.6 

When  MSOFE  [48]  begins,  one  of  the  many  files  it  reads  is  the  FLIGHT  file  provided  by 
PROFGEN  (See  Section  4.1).  Another  file  MSOFE  must  read  for  this  thesis  is  the  SV 
ephemeris  data  (provided  by  the  National  Geodetic  Office).  (This  thesis  is  the  first  at 
AFIT  to  use  true  SV  ephemeris  data  provided  by  the  National  Geodetic  Office).  The 
method  used  to  insert  the  true  ephemeris  data  into  MSOFE  is  by  modifying  a  FORTRAN 
routine  [57]  which  appends  true  ephemeris  data  to  the  end  of  the  binary  FLIGHT  file.  A 
copy  of  this  Fortran  routine  ("  ADDSV.for")  is  found  in  Appendix  E.  A  sample  of  the  true 
SV  ephemeris  data,  obtained  by  the  National  Geodetic  Office  is  also  found  in  Appendix 
E. 

The  SV  ephemeris  used  in  the  LSM  simulation  was  selected  based  upon  the  best  four  SV 
available  for  a  random  day  (i.e.,  21  May  94  used  in  this  thesis)  by  the  software  SEM3.6 
[18].  The  task  to  choose  which  4  SVs  to  use  was  determined  as  follows: 

•  Random  day  (21  May  94  )  selected. 

(GPS  week  749,  Day  of  year  141). 

•  GPS  Almanac  data  file  "05 1994A.AL3"  (See  Appendix  E)  was  obtained 
(downloaded)  via  the  Coast  Guard  Bulletin  Board  Services  (CGBBS). 

•  L AT/LONG/ ALT  along  the  Tanker  flight  profile  was  noted  and  entered 
in  SEM  3.6. 


Best  4  S  V  based  on  PDOP  algorithm 
5°  mask  angle 
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•  Scenario  duration:  2  hours  [Begin  time  04:00  UTC 
(Coordinated  Universal  Time  or  08:00  Eastern  time)] 

SEM  3.6  uses  the  identical  position  dilution  of  precision  (PDOP)  algorithm  that  the 
Collins  Receiver  3-A  P-Code  GPS  receiver  [13]  utilizes.  The  best  four  satellite  vehicles 
(SV)  are  then  numerically  displayed  to  the  user.  The  5°  mask  angle  was  chosen  based  on 
the  author's  knowledge  as  a  realistic  number  utilized  by  most  GPS  receivers.  "Mask 
angle"  is  defined  in  this  thesis  as  the  angle  made  from  the  surface  of  the  GPS  receiver 
antenna.  Thus  all  satellites  in  view  above  a  5°  angle  made  from  the  GPS  antenna  surface, 
on  21  May  94,  between  the  hours  of  04:00  -  06:00  GMT  are  "fair-game"  to  use  in  the 
selection  of  "Best-4"  SVs  based  on  a  PDOP  algorithm.  "GMT"  is  defined  as  Greenwich 
mean  time  (same  as  Zulu  or  Universal  time).  Plots  obtained  from  SEM3.6  (GDOP, 
PDOP,  VDOP,  TDOP,  SV  Bearing/Elevation,  SV  Rise/Set,  Number  of  Visible  SVs  and 
Elevation/Time)  are  shown  in  Appendix  F. 

The  Dilution  of  Precisions  (DOPs)  in  Appendix  F  are  defined  as: 

GDOP  =  ( a2  +  Oy  +  a:  +  a2  )1/2  =  Geometric  Dilution  of  Precision 
PDOP  -  (a2  +  a2  +  a2  )1/2  =  Position  Dilution  of  Precision 
HDOP  =  (a2  +  a2  )1/2  =  Horizontal  Dilution  of  Precision 

VDOP  =  (a2  )1/2  =  Vertical  Dilution  of  Precision 
TDOP  =  (a2  )1/2  =  Time  Dilution  of  Precision 

•  The  SV  Bearing/Elevation:  Displays  the  position  of  SVs  that  are  visible  at  the 
selected  location  during  the  specified  time  window  for  the  date  and  mask  angle 
chosen. 

•  SV  Rise/Set:  Shows  the  rise  times,  set  times,  visibility  periods,  and  selection  (in  use) 
periods  for  each  SV  as  functions  of  time. 

•  Number  of  Visible  SVs:  Shows  as  a  function  of  time,  the  total  number  of  SVs  visible 
during  the  time  period  under  evaluation. 

•  Elevation/Time:  Displays  the  elevation  of  visible  SVs  as  a  function  of  time. 

For  further  information  on  SEM  3.6  software  and  SEM  3.6  output  plot  format,  see  [18]. 
For  further  information  on  DOPs,  see  [17]. 

The  SVs  chosen  and  entered  in  the  LSM  simulation  are  summarized  in 
Table  4-5: 
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Time  (UTC) 

Satellite  Vehicle 
Number 

PROFGEN  elapsed 
time  (seconds) 

PDOP 

04:15 

6,16,17,28 

0 

2.25 

04:55 

6,12,17,28 

2400 

2.29 

05:10 

1,12,17,28 

3300 

2.47 

Table  4-5.  Satellite  Vehicles  Chosen  During  Simulation 


Note  that  changing  of  SVs  can  become  quite  evident  in  enhancing  or  degrading 
one's  overall  navigation  solution  (by  increasing  or  decreasing  DOP).  Sections  4.5  -  4.8 
clearly  show  evidence  of  degradation  when  a  SV  change  occurs  at  T  =  2400  seconds. 


4.3  Development  of  The  Three  INS's 

It  was  an  objective  in  this  thesis  to  compare  a  medium  accuracy  INS  (0.4  nm/hr)  with  a 
lower  quality  inertial  navigation  systems  in  the  range  of  2  to  4  nm/hr.  The  NRS  model 
(obtained  from  [45,50]  drifts  at  an  accuracy  of  approximately  0.4  nm/hr  (CEP).  This  drift 
rate  is  considered  medium  accuracy  for  this  thesis,  thus  two  lower  accuracy  INS  models 
were  created  (2.0  nm/hr  CEP  and  4.0  nm/hr  CEP).  The  choice  of  error  state  modification 
to  the  existing  NRS  INS  model  [45,50]  was  based  in  part  from  personal  conversations 
with  [31].  Modifications  to  the  system  covariance  matrices  for  each  respective  INS  are 
shown  in  Table  4-6.  Only  random  constant  shaping  filter  states  were  changed  (not  the 
1st  order  Gauss-Markov  drift  states,  etc.).  To  avoid  "CEP"  confusion,  CEP  will  be 
defined  next,  directly  from  [10]. 

In  the  case  of  two  jointly  Gaussian  variables  x  and  y,  where  ox=Gy=  a,  and 
p  =  0,  the  normal  distribution 


/x,y  (*.)0 


-1 

(x-x)2  2p(x-x)(y-y)  ( (y-y)2 

1 

_ ^  2(1— p2) 

a2  axay  a2 

2jtaxayV1-p2 


(4.1) 


is  called  circular  [10]. 


4-7 


LSM 

LSM  State  Name 

0.4  nm/hr 

2.0  nm/hr 

4.0 

SNU 

NRS 

Units 

State 

Value 

Value 

nm/hr 

84-1 

State 

Number 

Value 

State 

23 

X-component  of  gyro 

drift  rate  repeatability 

1.78E-16 

1 .6  IE- 14 

5.46E-14 

30 

25 

(°/hr)/hr 

24 

BBSS!! 

2. 1  IE- 16 

1.96E-14 

6.35E-14 

31 

26 

(°/hr)/hr 

25 

Z-component  of  gyro 
drift  rate  repeatability 

2.11E-16 

1.96E-14 

6.35E-14 

32 

27 

(7hr)/hr 

26 

X-component  of  gyro 
scale  factor  error 

2.49E-1 1 

9.09E-10 

3.00E-9 

33 

28 

% 

27 

Y -component  of  gyro 
scale  factor  error 

2.50E-1 1 

9.00E-10 

3.00E-9 

34 

29 

% 

28 

Z-component  of  gyro 
scale  factor  error 

2.49E- 1 1 

9.09E-10 

3.01E-9 

35 

30 

% 

29 

X-component  of 
accelerometer  bias 
repeatability 

2.31E-7 

2.17E-5 

7.26E-5 

48 

31 

M-g 

30 

Y -component  of 
accelerometer  bias 
repeatability 

2.31E-7 

2.17E-5 

7.26E-5 

49 

32 

M-g 

31 

Z-component  of 
accelerometer  bias 
repeatability 

2.30E-7 

2.16E-5 

7.25E-9 

50 

33 

M^g 

32 

X-component  of 
accelerometer  and 
velocity  quantizer  scale 
factor  error 

1.44E-8 

5.18E-7 

1.73E-6 

51 

34 

% 

33 

Y-component  of 
accelerometer  and 
velocity  quantizer  scale 
factor  error 

1.44E-8 

5.14E-7 

1.73E-6 

52 

35 

% 

34 

Z-component  of 
accelerometer  and 
velocity  quantizer  scale 
factor  error 

1.44E-8 

5.14E-7 

1.73E-6 

53 

36 

% 

35 

X-component  of 

accelerometer  and 
velocity  quantizer  scale 
factor  asymmetry 

6.25E-10 

7.50E-9 

2.50E-8 

54 

37  j 

% 

36 

Y-component  of 

accelerometer  and 
velocity  quantizer  scale 
factor  asymmetry 

6.25E-10 

7.55E-9 

2.50E-8 

55 

38 

% 

37 

Z-component  of 
accelerometer  and 
velocity  quantizer  scale 
factor  asymmetry 

6.25E-10 

7.50E-9 

2.50E-8 

56 

39 

% 

38 

X  accelerometer 
misalignment  about  Z- 
axis 

8.74E-15 

3.14E-13 

3.17E-12 

66 

40 

arc-sec 

39 

Y  accelerometer 
misalignment  about  Z- 
axis 

3.76E-10 

9.46E-13 

1.40E-7 

67 

41 

arc-sec 

40 

Z  accelerometer 
misalignment  about  Y- 
axis 

3.76E-10 

9.40E-13 

1.40E-7 

68 

42 

arc-sec 

41 

X  accelerometer 
misalignment  about  Y- 
axis 

3.76E-10 

9.40E-13 

1.40E-7 

69 

43 

arc-sec 

4-8 


Table  4-6.  MSOFE  Initial  System  Covariance  Parameters  for  0.4/2. 0/4.0  nm/hr  INS 
where 


x  =  £[x] 
y  =  E[y] 

a2=£[(x-x)2]  (4.2) 

a2=Z?[(y-y)2] 

p  =  £[(x-x)(y-y)/axay] 

(Here:  p  =  the  correlation  coefficient  of  x  and  y,  ox  is  the  standard  deviation  of  x  and  oy 

is  the  standard  deviation  of  y).  Then  the  probability  density  for  radial  distance  r  from  the 
point  (x,y)  corresponding  to  f(x,y)  is 

f(r)  =  (r/a2)e~rl/2a2  (4.3) 


where 


r2  =  (jc  — x)2  +(y  — y)2  (4.4) 

r  is  called  the  radial  error  [10].  The  equal  probability  ellipses  are  circles  in  this  case, 
since  or  =  av  =  a.  The  probability  p'  that  a  point  (x,y)  taken  at  random  will  fall  within  a 

x  y 

circular  ring,  whose  center  is  at  (x,y)  and  whose  inner  and  outer  radii  are  r  and  r  +  Ar,  is 

p'  =  f(r)Ar  (4.5) 

The  probability  p  that  a  point  (x,y)  taken  at  random  will  fall  within  the  circle 

(x-x)2  +(y-y)2  =  (cg)2  (4.6) 


for  a  particular  value  of  c  is 


P  = 


-„-V2 


(4.7) 
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When  c  =  1. 1774,  p  =  0.5,  and  the  circle  of  radius  1 . 1774a  is  called  the  50  percent 
probability  circle  [10],  The  radius  of  this  circle  is  called  the  circular  probable  error  or 
CEP  and 


CEP  =  1.1774a  (4.8) 

Approximations  are  often  made  for  the  CEP  calculation  when  the  actual  error 
distribution  has  elliptical  rather  than  circular  equal-probability  loci.  The  approximation 
used  in  this  thesis  is: 


Estimate(CEP)  =  0.5881a*  -t-a^,]  (4.9) 

The  CEP  curves  for  one  hour  of  operation  of  the  three  inertial  navigation  systems  created 
in  this  thesis  are  shown  in  Figure  4  -  1 . 
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4.4  Filter  Tuning  Example 


This  section  describes  the  author's  approach  to  tuning  the  filter  error  states.  Obviously,  to 
land  an  aircraft  safely,  a  3-dimensional  representation  of  position  is  of  utmost  importance, 
particularly  with  regard  to  the  vertical  axis.  Also  of  importance  to  the  author  was  to 
estimate  the  GPS  clock  bias  error  adequately  since  this  is  the  largest  contributor  to  the  P- 
Code  GPS  receiver  error.  Again,  looking  at  the  "big-picture",  the  truth  model  used  in  the 
simulation  is  composed  of  69  error  states  which  attempt  to  represent  the  actual  error 
outputs  of  the  real-world  avionics  on  board  our  simulated  aircraft.  Theoretically,  69  error 
states  is  a  minuscule  amount  when  actually  an  infinite  amount  of  error  states  would  be 
needed  to  represent  the  real-world  avionics;  nevertheless,  the  69  error  states  are  assumed 
to  provide  an  adequate  truth  source. 

Our  truth  model  (numbering  69  error  states)  requires  a  large  amount  of  computer 
processing,  and  a  filter  based  upon  this  truth  model  would  not  be  feasible  at  the  present 
time  to  use  onboard  civil  aviation  aircraft.  Our  objective  in  this  thesis  is  to  minimize 
computer  processor  time  (through  number  of  error  states)  and  costs,  thus  we  try  to  use  a 
subset  of  the  truth  model  error  states  as  the  filter  model  error  states.  If  we  are  successful, 
we  will  be  able  to  choose  a  truth  model  subset  (in  this  thesis  called  the  filter  model, 
which  is  composed  of  13  states)  that  adequately  models  the  truth  model  (remembering 
that  the  truth  model  represents  the  real  world).  We  will  add  a  specific  amount  of  white 
noise  to  the  differential  equations  for  our  individual  filter  model  error  states  to 
compensate  for  truth  model  error  states  we  were  unable  to  model  specifically  and 
distinctly  using  our  "less-complex"  filter  model.  Performing  this  procedure  is  considered 
tuning  the  dynamics  noise  strengths  or  "Q"  values  of  the  filter  model.  (Another  parameter 
tuned  is  the  "R"  value  of  the  measurement  equations,  i.e.,  the  measurement  noise 
covariance  matrix.  The  "R"  tuning  is  not  described  in  this  example,  though  it  follows 
similar  procedures.  For  more  information  see  [40,41,42]).  Recall  "Q"  and  "R"  from 
earlier  model  equations  of  Chapter  3. 

Let's  look  at  several  states  of  Case  I  (which  utilizes  a  0.4  nm/hr  INS)  that  were 
tuned  using  the  legend  in  Figure  4-2. 
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Plot  Legend: 

true  error  (sample  mean  error  ±  ctme) 
filter-predicted  error  (0  ±  o^iter) 
—  sample  mean  error 

Figure  4-2.  Plot  Legend 


The  plots  contained  in  Appendix  G  -  Appendix  U  contain  five  traces.  The  innermost 

trace  (a  solid  line - )  on  each  data  plot  is  the  sample  mean  error  time  history  for  the 

applicable  error  state  and  is  defined  by  [58,40]: 


1  N 


(4.10) 


where  Xj(tt)  is  the  filter-computed  estimate  of  a  given  navigation  variable  and  xtruej  (tt ) 
is  the  truth  model  value  of  the  same  variable,  at  time  tt ,  for  sample  run  j.  N  is  the  number 
of  Monte  Carlo  runs  in  the  simulation  (15  in  this  thesis)  [45]. 


In  addition  to  the  center  trace,  two  more  pairs  of  lines  are  plotted.  The  first  pair 

(represented  by  "dots"  . in  the  plots)  is  the  Mean± Sigma.  The  Mean± Sigma  is 

symmetrically  displaced  about  the  mean  error,  Me(ti).  The  Mean  ±  Sigma  is  the  sum  of 
the  previously  defined  mean,  Me(f(),  and  the  actual  filter  standard  deviation  ±^jPe(ti) , 
where  Pgitf)  is  the  true  error  sample  variance  at  time  t{ .  The  true  standard  deviation  is 
calculated  from  the  following  equation  [45,40]: 

<W<i)  =  =  tfw  (4.1 1) 

where  N  is  the  number  of  runs  in  the  Monte  Carlo  simulation  (15  in  this  thesis),  and 
Mg(ti)  is  the  square  of  the  mean  of  the  variable  at  each  time  of  interest. 


The  last  pair  of  traces  (represented  by  "dashes" - )  is  the  filter-computed  ±G 

filter  values  for  the  same  variables  of  interest  and  are  symmetrically  displaced  about  zero 
because  the  filter  "believes"  that  it  is  producing  zero  mean  errors  [45].  These  quantities 
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are  propagated  and  updated  in  the  MSOFE  [48]  software  using  the  covariance 
propagation  equation  shown  in  Chapter  2.  Thus,  this  last  pair  of  traces  represents  the 
filter's  estimate  of  the  size  of  its  own  error. 


We  first  begin  this  example  with  tuning  the  three  tilt  error  states.  A  fruitful  tuning 
strategy  was  found  to  entail  tuning  the  platform  (tilts)  first,  then  tuning  the 
accelerometers.  Looking  at  Figure  G  -  1  in  Appendix  G,  we  see  conservative  tuning  for 
the  North,  West  and  Azimuth  Tilts  when  compared  with  Figure  G  -  2  in  Appendix  G  (the 
true  error,  mean  error  ±  alrue,  is  well  within  filter-computed,  0  ±  cfilter,  error).  Figure  G  - 
2  has  been  tuned  by  decreasing  the  Q  values  as  shown  in  Table  4-7. 


LSM  Filter 

State 

Number 

LSM  State 

Name 

Q- value 
(Before) 
(ft2/s4) 

Q- value 
(After) 
(ft2/s4) 

SNU-84 

State 

Number 

NRS  State 
Number 

4 

N-Tilt 

9.52E-13 

8.56E-13 

4 

4 

5 

W-Tilt 

9.52E-13 

8.56E-13 

5 

5 

6 

Az-Tilt 

1.62E-11 

1.52E-11 

6 

6 

Table  4-7.  Q-Values  for  N,  W  and  Az  Tilts  (Before  and  After  Tuning) 


The  North,  West  and  Vertical  velocity  error  states  were  tuned  conservatively  at 
the  start  as  shown  in  Figure  G  -  3.  Figure  G  -  4  shows  final  tuning  of  North,  West  and 
Vertical  velocities.  Note  that  the  vertical  velocity  plot  in  Figure  G  -  4  may  not  be 
conservative  enough:  true  mean  error  ±  atrue  does  not  always  stay  within  0  ±  Oyjter 
Table  4-8  shows  the  old  and  new  filter  tuning  Q  values  of  Figure  G  -  3  and  Figure  G  -  4. 


LSM  Filter 
State 
Number 

LSM  State 
Name 

Q- value 
(Before) 
(ft2/s4) 

Q- value 
(After) 
(ft2/s4) 

SNU-84 

State 

Number 

NRS  State 
Number 

7 

BIS 

1.54E-3 

5.15E-6 

7 

7 

8 

1.54E-3 

5.15E-6 

8 

8 

9 

7.20E-3 

7.72E-3 

9 

9 

Table  4-8.  Q-Values  for  N,  W  and  Vertical  Velocity  (Before  and  After  Tuning) 


Lastly,  the  remainder  of  the  filter  states  were  tuned.  The  results  of  these  filters 
tunings  are  used  for  Case  I.  If  the  reader  is  interested,  final  tuning  plots  for  all  filter  states 
of  this  example  can  be  found  by  looking  in  Appendix  H,  Case  I  tuning  plots. 


4-13 


4.5  Radar  Altimeter/Baro/GPS  Aiding  of  the  0.4  nm/hr,  2.0  nmlhr  and  4.0  nm/hr 

INS 

This  section  discusses  the  results  of  Cases  I- VI.  Case  I- VI  plots  are  located  in  Appendix 
H  -  Appendix  M.  Cases  I,  in  and  V  do  not  use  radar  altimeter  aiding,  whereas.  Cases  II, 
IV  and  VI  do. 

4.5.1  Case  I  and  Case  II  -  0.4  nm/hr  INS 

•  Case  I  (Figures  H  -  1  -  H  -  6  in  Appendix  H) 

The  Tanker  flight  profile  as  described  in  Section  4. 1  was  flown  with  no  outages 
(at  all  times)  of  4  SV  GPS  coverage.  Throughout  the  flight,  the  INS  errors  were  bounded 
by  GPS  aiding.  Filter-predicted  error  of  longitude  during  the  last  1 1  minutes  of  flight 
slightly  underestimate  the  true  error  (Figure  H-l  and  Figure  H-3  being  of  particular 
interest  with  respect  to  performance).  This  slight  underestimation  (2  feet)  can  easily  be 
corrected  by  adjusting  the  'Q'  value  for  the  respective  horizontal  error  state  following  the 
tuning  example  described  in  Section  4.4.  The  primary  goal  of  this  thesis  was  to  show 
filter  performance  trends;  'Q'  value  filter  tuning  as  shown  for  all  cases  is  reasonable, 
though  not  in  an  ultimate  final  form  for  flight  test  purposes  due  to  time  constraints.  Note 
again  that,  for  all  cases,  changing  of  SVs  can  become  quite  evident  in  enhancing  or 
degrading  one's  overall  navigation  solution  (by  increasing  or  decreasing  DOP).  Sections 
4.5  -  4.8  clearly  show  evidence  of  degradation  when  a  SV  change  occurs  at  T  =  2400 
seconds. 

Note  that  conservative  filter  performance  of  Table  4  -  9  is  due  to  conservative 
Q  filter vs-  ® true  R'  tuning  of  the  satellite  vehicle  measurements.  The  a  fllter  and  atrue 

values  chosen  for  satellite  vehicle  'R'  measurements  were  8.66  ft  and  1.414  ft, 

respectively;  See  Table  B  -  1 1  in  Appendix  B  (discussed  in  Section  I  for  Case  II).  The 
choices  of  and  otrue  dictate  a  "flooring"  value  seen  in  the  true  and  filter 

performance.  The  'R'  tuning  values  stated  above  are  used  for  all  cases  in  Section  4.5  - 
4.8;  thus  in  all  cases  investigated  in  this  thesis,  the  filter  performance  will  show  overly 
conservative  results  (evident  in  Tables  4-9, 4  - 11,4  -  13,  and  4-21). 

.  Case  13  (Figures  I-l-I-6in  Appendix  I) 

This  is  the  same  as  Case  I  except  at  approximately  3000  ft  (AGL),  the  radar 
altimeter  begins  to  provide  measurements  to  the  Kalman  filter.  The  accuracy  of  the  radar 
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altimeter  measurements  increases  (i.e.  "R”  value  decreases)  as  altitude  above  the  ground 
decreases  (see  Section  3.4.2,  radar  altimeter  description). 

The  following  error  states  showed  improvement  using  the  radar  altimeter 
measurements: 


Aircraft  altitude  (See  Figure  1-2  and  1-3,  compared  to  Figures  H  -  2  and 
H  -  3) 

Longitude  error  (See  Figure  I  -  1,  compared  to  H  -  1) 

GPS  user  clock  bias  (See  Figure  1-6,  compared  to  H  -  6) 


Table  4-9  summarizes  Case  I  and  Case  II  landing  system  performance  at  respective 
decision  heights. 


Position  Error 
State 

Decision 
Height  (feet) 

Case  I  True 
Error  (feet) 

Case  I  Filter 
Error  (feet) 

Case  II  True 
Error  (feet) 

Case  II  Filter 
Error  (feet) 

Slatitude 

200 

8.81 

9.85 

8.72 

9.69 

Slatitude 

100 

9.40 

9.80 

9.13 

9.69 

Slatitude 

50 

8.61 

9.80 

8.35 

9.69 

Slongitude 

200 

9.14 

10.68 

8.83 

10.58 

Slongitude 

100 

9.77 

10.68 

9.65 

10.58 

50 

8.83 

10.69 

8.70 

10.58 

Saltitude 

200 

14.2 

15.1 

3.76 

7.74 

Saltitude 

100 

15.65 

15.1 

2.43 

7.57 

8altitude 

50 

15.3 

15.1 

2.64 

7.52 

Table  4-9.  Case  I  and  Case  II  la  Latitude,  Longitude  and  Altitude  Errors 


Clearly,  as  indicated  in  the  last  3  rows  of  Table  4-9,  the  radar  altimeter  reduces  the 
altitude  error.  A  zoomed-in  plot  of  horizontal  errors  latitude,  longitude  and  vertical  errors 
using  the  radar  altimeter  for  Case  II  is  found  in  Figure  I  -  3  of  Appendix  I.  "Flooring" 
occurs  in  the  altitude  plot  of  Figure  I  -  3  as  the  aircraft  gets  lower  in  altitude;  the  filter 
becomes  more  overly  conservative.  Re-tuning  to  achieve  better  tuning  for  all  time, 
especially  near  decision  points,  would  be  useful  to  accomplish.  Note  the  increased 
accuracy  when  the  radar  altimeter  measurements  become  available  at  T  >  3641  seconds. 
The  radar  altimeter  (Case  II)  reduced  the  aircraft  altitude  error  during  the  final  approach 
by  approximately  80%  (true)  and  49%  (filter-predicted). 

Table  4-10  summarizes  the  Case  I  and  Case  II  performance  in  terms  of  what 
Category  of  Precision  Approach  is  met  based  on  Horizontal  and  Vertical  errors. 
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Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 

(Based  on  Filter 
Predicted  Error) 

I 

Cat  I 

Cat  I 

None 

None 

II 

Cat  I 

Cat  I 

Cat  II 

None 

Table  4-10.  Precision  Landing  Category  Predicted  To  Achieve 


The  filter-predicted  error  in  the  vertical  channel  with  radar  altimeter  aiding  exceeds  the 
Category  I  landing  specification  by  approximately  1  foot,  but  the  true  error  shows  that  the 
system  with  the  radar  altimeter  could  possibly  meet  Category  II  landings  vertically,  but 
not  Category  II  horizontally.  For  Case  I  (Baro-al  time  ter  and  P-Code  GPS  aiding)  a 
precision  approach  is  impossible  due  to  the  large  vertical  errors.  Additional  tuning 
('Rfiiter')  °f  Case  II  filter  is  required  to  meet  a  precision  approach,  but  it  was  not 
accomplished  in  this  thesis  due  to  time  constraints. 


4.5.2  Case  III  and  Case  IV  -  2  nm/hr  INS 
.  Case  III  (Figures  J  -  1  -  J  -  6  in  Appendix  J) 

The  Tanker  profile  as  described  in  Section  4.1  was  flown  with  no  outages  of  4  SV  GPS 
(coverage  at  all  times).  Throughout  the  flight,  the  2.0  nm/hr  INS  errors  were  bounded  by 
GPS  aiding.  Filter-predicted  error  of  longitude  for  the  first  1 1  minutes  of  flight 
underestimates  the  true  error.  Adjustment  of  the  'Q'  value  for  the  respective  horizontal 
error  state  should  be  implemented  as  described  in  Section  4.4,  though  not  done  due  to 
time  constraints. 

•  Case  IV  (Figures  K  -  1  -  K  -  6  in  Appendix  K) 

This  is  the  same  as  Case  III  except  at  approximately  3000  ft  (AGL),  the  radar 
altimeter  begins  to  provide  measurements  to  the  Kalman  filter.  The  accuracy  of  the  radar 
altimeter  measurements  increases  as  altitude  decreases.  Using  the  radar  altimeter  (Case 
IV)  reduced  the  aircraft  altitude  error  during  final  approach  by  approximately  84%  (true) 
and  47%  (filter-predicted). 

Table  4-11  summarizes  Case  in  and  Case  IV  landing  system  performance  at 
respective  decision  heights.  Table  4-12  summarizes  the  Case  HI  and  Case  IV 
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performance  in  terms  of  what  Category  of  Precision  Approach  is  met  based  on 
Horizontal  and  Vertical  errors. 


In  Case  III,  only  the  horizontal  components  of  position  will  meet  Category  I 
requirements.  Case  IV  results  show  that  the  current  filter  overestimates  the  true  error  by 
244%  at  Decision  Height  (DH)  of  200,  161%  at  DH  of  100  and  345%  at  DH  of  50. 
Additional  filter  tuning  specifically  during  the  landing  phase  of  the  flight  profile  is 
needed  to  meet  a  CAT  I  precision  approach.  Additional  tuning  was  not  performed  due  to 
time  constraints  of  this  project.  Conservative  tuning  of  the  filter  is  discussed  in  Chapter 
5. 


Position  Error 
State 

Decision 
Height  (feet) 

Case  III  True 
Error  (feet) 

Case  III  Filter 
Error  (feet) 

Case  IV  True 
Error  (feet) 

Case  IV  Filter 
Error  (feet) 

Slatitude 

200 

8.49 

9.90 

8.47 

9.75 

Slatitude 

100 

9.01 

9.90 

9.46 

9.74 

Slatitude 

50 

8.67 

9.90 

8.95 

9.74 

Slongitude 

200 

9.16 

10.75 

9.08 

10.60 

Slongitude 

100 

9.81 

10.76 

9.64 

10.65 

Slongitude 

50 

8.79 

10.76 

9.39 

10.65 

Saltitude 

200 

15.1 

15.8 

2.43 

8.37 

Saltitude 

100 

15.93 

15.8 

3.14 

8.21 

Saltitude 

50 

15.3 

15.8 

1.83 

8.16 

Table  4-11.  Case  III  and  Case  IV  la  Latitude,  Longitude  and  Altitude  Errors 


Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 
(Based  on  Filter 
Predicted  Error) 

III 

Cat  I 

Cat  I 

None 

None 

IV 

Cat  I 

Cat  I 

Cat  I 

None 

Table  4  -12.  Precision  Landing  Category  Predicted  To  Achieve 


4.5.3  Case  V  and  Case  VI  -  4  nmfhr  INS 

•  Case  V  (Figures  L  -  1  -  L  -  6  in  Appendix  L) 

•  Case  VI  (Figures  M  -  1  -M  -  6  in  Appendix  M) 

The  4.0  nm/hr  INS  was  bounded  by  GPS  for  Case  V  and  by  the  combination  of 
GPS/Radar  altimeter  measurement  for  Case  VI.  Again,  filter-predicted  error  of  Longitude 
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for  the  last  1 1  minutes  of  flight  may  underestimate  the  true  Longitude  error.  Adjustment 
of  the  Q  value  for  the  respective  horizontal  error  state  should  be  implemented  as 
described  in  Section  4.4.  Using  the  radar  altimeter  (Case  VI)  reduced  the  aircraft  altitude 
error  during  final  approach  by  approximately  85%  (true)  and  45%  (filter  predicted). 


Table  4-13  summarizes  Case  V  and  Case  VI  landing  system  performance  at  respective 
decision  heights. 


Position  Error 

Decision 

Case  I  True 

Case  I  Filter 

Case  II  True 

Case  II  True 

State 

Height  (feet) 

Error  (feet) 

Error  (feet) 

Error  (feet) 

Error  (feet) 

Slatitude 

200 

8.63 

9.95 

8.40 

9.79 

Slatitude 

100 

9.43 

9.95 

9.16 

9.78 

Slatitude 

50 

8.73 

9.95 

8.85 

9.80 

Slongitude 

200 

9.12 

10.8 

8.77 

10.70 

Slongitude 

100 

9.91 

10.8 

9.93 

10.71 

Slongitude 

50 

8.72 

10.8 

9.44 

10.71 

Saltitude 

200 

14.96 

15.75 

2.39 

8.23 

Saltitude 

100 

16.5 

15.75 

2.82 

8.10 

Saltitude 

50 

15.3 

15.75 

1.75 

8.00 

Table  4-13.  Case  V  and  Case  VI  la  Latitude,  Longitude  and  Altitude  Errors 


Table  4-14  summarizes  the  Case  V  and  Case  VI  performance  in  terms  of  what  Category 
of  Precision  Approach  is  met  based  on  Horizontal  and  Vertical  errors. 


Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 
(Based  on  Filter 
Predicted  Error) 

V 

Cat  I 

Cat  I 

None 

None 

VI 

Cat  I 

Cat  I 

Cat  I 

None 

Table  4-14.  Precision  Landing  Category  Predicted  To  Achieve 


Additional  filter  tuning  specifically  during  the  landing  phase  of  the  flight  profile  is 
needed  to  meet  a  CAT  I  precision  approach.  Due  to  time  constraints,  final  filter  tuning 
was  not  accomplished,  but  will  be  discussed  as  an  issue  in  Chapter  5. 
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4.6  Radar  Altimeter/Baro/GPS  Aiding  of  the  0.4  nm/hr,  2.0  nmfhr  and  4.0  nm/hr 
INS  using  an  Additional  Single  Pseudolite  during  the  Landing  Approach 

This  section  discusses  the  enhancement  when  an  additional  measurement  from  a  single 
pseudolite  is  positioned  along  the  final  approach  path.  Averaged  results  from  Cases  VII  to 
XII  will  be  compared  with  the  baseline  Cases  I,  in  and  V  of  Section  4.5.  "Averaged" 
values  were  permitted  to  be  used  based  on  the  fact  that  no  significant  performance 
difference  is  noted  between  each  respective  INS  when  there  were  no  GPS  measurement 
outages.  An  illustration  of  the  aircraft  final  approach  and  pseudolite  location  is  found  in 
Section  4.6.1. 

4.6.1  Cases  VII,  VIII  and  IX  -  Pseudolite  But  No  Radar  Altimeter 

Case  VII,  Case  VIII  and  Case  EX  follow  the  configuration  found  in  Table  4-15. 


Case  VII 

Case  Vm 

Case  IX 

Barometric 

Barometric 

Barometric 

Altimeter 

Altimeter 

.  l 

0.4  nm/hr 

2.0  nm/hr 

4.0  nm/hr 

CEP  INS 

CEP  INS 

CEP  INS 

P-Code 

P-Code 

P-Code 

GPS 

GPS 

GPS 

Pseudolite 

Pseudolite 

Pseudolite 

Table  4-15.  Case  Vn,  Case  Vffl  and  Case  IX  Configurations 

•  Case  VII  (Figures  N-l-N-6in  Appendix  N) 

•  Case  VHI  (Figures  O-  1  -  O  -  6  in  Appendix  O) 

•  Case  IX  (Figures  P  -  1  -P  -  6  in  Appendix  P) 

The  pseudolite  placement  was  chosen  to  minimize  VDOP  errors  during  a  portion 
of  the  aircraft  landing.  VDOP  decreased  from  1.7  to  1.3  using  the  single  pseudolite.  This 
decrease  in  VDOP  (23  %)  is  reflected  in  the  decrease  in  altitude  error  (21.3%)  seen  when 
comparing  the  averaged  results  of  Cases  I,  El,  V  vs.  Cases  VII  -  IX;  see  Table  4-16.  It 
should  be  noted  that  the  GPS  receiver  maintains  lock  on  the  best  four  SV  overhead,  while 
maintaining  lock  on  the  fifth  SV  (the  pseudolite  on  the  ground).  The  designer  must 
maintain  the  insight  on  importance  of  geometry  if  it  is  necessary  to  replace  one  of  the 
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existing  SVs  overhead,  with  a  ground  based  pseudolite  instead.  The  GPS  receiver  must 
simultaneously  maintain  an  adequate  (minimal)  HDOP  and  PDOP  at  all  times  during  the 
precision  approach.  Alternate  locations  of  pseudolites  and  their  impact  on  performance 
(due  to  geometry)  must  be  chosen  with  care. 

The  aircraft  is  assumed  to  be  following  a  precision  approach  pattern,  landing  at 
Wright-Patterson  AFB,  OH  as  found  in  [15].  A  pseudolite  is  placed  at  the  coordinates 
indicated  in  Figures  4-3  and  4-4.  At  T=3641,  the  pseudolite  measurement  is  available 
and  being  processed  once  every  second.  The  aircraft  has  continued  its  descent  from  4028 
feet  (MSL)  to  3799  feet  (MSL)  (where  elapsed  time  is  now,  T=3658  seconds)  and  makes 
a  +45  turn  to  align  itself  with  the  runway  (Heading  232°).  The  aircraft  continues  its 
descent  at  a  3°  glideslope,  now  aligned  with  the  runway  (Heading  232°)  at  a  velocity  of 
133  knots.  The  aircraft  lands  at  an  elapsed  time  of  T=3913  seconds  true  altitude  is  0  ft 
(AGL),  825  ft  (MSL). 

Looking  at  the  plots  of  Cases  VII,  VIE  and  IX  (in  Appendix  N,  Appendix  O  and 
Appendix  P),  it  is  clearly  evident  that  the  Latitude,  Longitude,  Altitude  and  GPS  clock 
bias  errors  are  significantly  reduced.  Seen  at  T  =  2400  and  T  =  3300  are  large  effects  of 
latitude,  longitude  and  clock  bias  error  due  to  satellite  vehicle  changes.  Table  4-16 
shows  the  average  true  error  reduction  by  using  the  single  pseudolite  when  compared 
when  no  pseudolite  is  used  (Case  I,  Case  HI  and  Case  V).  Note  again  that  Cases  I,  HI  and 
V  results  are  averaged  due  to  similar  results  between  INS  types  with  no  GPS  outages. 


Figure  4-3.  Top  View  of  Tanker  Precision  Approach 
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3912  3800  3658  3641 

Elapsed  Flight  Time  (seconds) 


Figure  4-4.  Side  View  of  Tanker  Precision  Approach 


Error  State 

Case  I,  III,  V  average 
true  error  (feet) 

Case  VII  -  IX  average 
true  error  (feet) 

%  Change  in  Error 

Latitude 

8.9 

4.6 

48.3% 

Longitude 

9.2 

3.9 

57.6% 

Altitude 

15.0 

11.8 

21.3% 

Clock  Bias 

9.0 

4.9 

45.6% 

Table  4-16.  Averaged  True  Error  Reduction  Using  a  Single  Pseudolite 


It  should  be  remembered  that  by  no  means  is  the  location  of  the  single  pseudolite 
"optimal".  No  real  criterion  was  used  when  selecting  the  pseudolite  location.  As  can  be 
seen  in  Figure  4  -  4,  by  the  time  the  aircraft  is  near  the  Category  I,  200-foot  decision 
height,  the  pseudolite  location  will  primarily  aid  in  the  horizontal  channel.  Pseudolite 
placement  nearer  to  the  decision  heights  may  aid  the  vertical  channel  by  an  additional 
12%  to  23%;  the  goal  is  to  bring  and  maintain  VDOP  from  1.3  to  1.0.  Readers  may  wish 
to  consult  [12]  for  further  reference  on  pseudolites.  Also,  it  was  assumed  no  multipath 
errors  exist  at  the  given  location  and  SV  clock  bias,  SV  position  errors,  tropospheric  and 
ionospheric  errors  could  be  ignored.  Thus,  for  an  aircraft  using  a  0.4,  2.0  or  4.0  nm/hr 
INS,  using  a  P-code  receiver  (no  differential  corrections  or  carrier-phase  measurements) 
and  the  assumed  location  (non-optimal)  of  a  single  pseudolite,  Category  I  precision 
approach  is  met  for  horizontal  errors,  but  no  category  of  precision  approach  is  possible  in 
the  vertical  axis  (again,  based  on  the  non-optimal  location  of  the  pseudolite). 
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Table  4-17  summarizes  the  Case  VII,  Case  VIII  and  Case  XI  performance  in  terms  of 
what  category  of  precision  approach  could  be  met  based  on  horizontal  and  vertical  errors. 


Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 
(Based  on  Filter 
Predicted  Error) 

VII-IX 

Cat  I 

Cat  I 

None 

None 

Table  4-17.  Precision  Landing  Category  Predicted  To  Achieve 


4.6.2  Cases  X,  XI  and  XII  -  Pseudolite  and  Radar  Altimeter 
Cases  X,  XI  and  XII  follow  the  configuration  found  in  Table  4-18. 


•  Case  X  (Figures  Q  -  1  -  Q  -  6  in  Appendix  Q) 

.  Case  XI  (Figures  R-  1  -  R  -  6  in  Appendix  R) 

•  Case  XII  (Figures  S  -  1  -S  -  6  in  Appendix  S) 

Cases  X,  XI  and  XII  are  identical  to  the  three  cases  of  the  previous  section  (Section  4.6.1) 

except  at  an  elapsed  flight  time  of  T=3800  (aircraft  on  final  approach  at  ALT  of  2139 
MSL  or  1314  AGL)  the  radar  altimeter  measurements  are  available.  The  combination  of 
the  pseudolite  and  radar  altimeter  give  performance  enhancements  summarized  in 
Table  4  -  19  to  Table  4-21.  A  large  improvement  can  clearly  be  seen  in  Latitude, 
Longitude  and  GPS  clock  bias  when  the  radar  altimeter  measurements  are  available  at  T  s 
3800  seconds  in  Appendices  Q,  R  and  S. 


Case  Xn 

Barometric 

Barometric 

Barometric 

Altimeter 

Altimeter 

Altimeter 

0.4  nm/hr 

2.0  nm/hr 

4.0  nm/hr 

CEP  INS 

CEP  INS 

CEP  INS 

P-Code 

P-Code 

GPS 

GPS 

Radar 

Radar 

Radar 

Altimeter 

Altimeter 

Altimeter 

and 

and 

and 

Pseudolite 

Pseudolite 

Pseudolite 

No  GPS 

No  GPS 

No  GPS 

Outage 

Outage 

Outage 

Table  4  -  18.  Case  X,  Case  XI  and  Case  XII  Configurations 
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Error  State 

Case  I,  III,  V  average 
true  error  (feet) 

Case  X,  XI,  XII  average 
true  error  (feet) 

%  Change  in  Error 

Latitude 

8.9 

4.4 

50.6% 

Longitude 

9.2 

3.3 

64.1% 

Altitude 

15.0 

2.4 

84.0% 

Clock  Bias 

9.0 

4 

55.0% 

Table  4  -  19.  Averaged  True  Error  Reduction,  Case  I-H-III  vs.  Case  X-XI-XII 


Error  State 

Case  II,  IV,  VI  average 
true  error  (feet) 

Case  X  -  XII  average 
true  error  (feet) 

%  Change  in  Error 

Latitude 

8.8 

4.4 

50.0% 

Longitude 

9.2 

3.3 

64.1 

Altitude 

2.6 

2.4 

7.7% 

Clock  Bias 

7.5 

4 

46.7% 

Table  4  -  20.  Averaged  True  Error  Reduction,  Case  II-IV-VI  vs.  Case  X-XI-XII 


Error  State 

Case  VII,  VIII,  IX 
average  true  error  (feet) 

Case  X  -  XII  average 
true  error  (feet) 

%  Change  in  Error 

Latitude 

4.6 

4.4 

4.4% 

Longitude 

3.9 

3.3 

15.4% 

Altitude 

11.8 

2.4 

79.7% 

Clock  Bias 

4.9 

4 

18.4% 

Table  4-21.  Averaged  True  Error  Reduction,  Case  VII-Vin-IX  vs.  Case  X-XI-XII 


Table  4-22  summarizes  the  Case  X,  Case  XI  and  Case  XII  performance  in  terms  of  what 
category  of  precision  approach  could  be  met  based  on  horizontal  and  vertical  errors. 


Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 
(Based  on  Filter 
Predicted  Error) 

X  -  XII 

Cat  II 

Cat  I 

Cat  II 

None 

Table  4  -22.  Precision  Landing  Category  Predicted  To  Achieve 


4.7  Example  of  0.4  nm/hr,  2.0  nm/hr  and  4.0  nm/hr  INS  Drift  during  a  Loss  of 
GPS  Measurement  for  1/8  Schuler  Period. 

•  Cases  Xm,  XTV  and  XV  (Figures  T  -  1  -T  -  13  in  Appendix  T) 
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Earlier  cases  (Cases  I-XII,  Sections  4.5  and  4.6)  showed  very  little  performance 
differences  when  comparing  say,  a  0.4  nm/hr  INS  and  a  4.0  nm/hr  INS.  Performance 
differences  were  noted  as  slight  because  no  degradation  to  the  baro-altimeter,  P-code  GPS 
or  radar  altimeter  measurements  occurred.  This  section  examines  the  possibility  of  loss 
of  GPS  for  approximately  650  seconds  (1/8  Schuler  period).  This  experiment  was  chosen 
based  upon  noting  that,  in  prior  Cases  I  -  XU,  the  three  INS's  do  not  seem  to  yield  very 
different  results  with  continuous  GPS  aiding.  In  the  real  world,  loss  of  GPS 
measurements  will  unfortunately  be  the  case  and  therefore  loss  of  GPS  measurements 
was  investigated  to  view  the  performance  of  the  true  and  filter  parameters  under  a  short 
time  duration  of  GPS  outage.  The  1/8  Schuler  period  was  chosen  so  it  could  be  possible 
to  multiply  the  error  state  degradations  by  a  factor  of  8,  if  one  wishes  to  simulate  a  full 
84-minute  Schuler  period  (also,  the  flight  profile  totals  only  65  minutes). 

Basically,  at  an  elapsed  time  of  T=2000  seconds  of  the  Tanker  profile,  loss  of  4 
SV  measurements  occurs.  Loss  of  4  SV  simulates  a  possible  flight  scenario  of  terrorist 
jamming,  an  incapacitated  GPS  receiver,  etc.  The  true  and  filter  errors  are  seen  to  grow 
larger  over  time,  until  the  GPS  measurements  are  re-acquired  at  T=2630  sec.  Once  GPS 
measurements  are  re-acquired,  estimation  of  the  INS  error  occurs.  Note  the  larger  true 
position  errors  of  the  4.0  nm/hr  system  (Figure  T-l,  bottom  plot)  vs.  the  smaller  true 
position  errors  of  the  0.4  nm/hr  system  (Figure  T  -  1,  top  plot);  be  aware  of  the 
dramatically  different  scales  for  the  three  plots.  Also  note  that  the  altitude  error  (Figure 
T  -  3)  grows,  but  becomes  bounded  by  baro  altimeter  aiding. 

The  point  of  this  section  is  to  show  that,  if  a  low  accuracy  INS  is  used  in  a 
GPS/INS  Kalman  filter  integration  for  a  landing  system  implementation,  one  must  realize 
the  consequences  if  the  GPS  measurement  data  is  unavailable,  especially  during  final 
approach.  The  results  of  this  section  are  directly  comparable  to  Cases  I,  ID  and  V.  Loss 
of  GPS  measurements,  but  with  a  system  supplemented  by  a  radar  altimeter  is  explored 
further  in  Section  4.8. 

The  plots  of  Appendix  T  also  show  that  the  4.0  nm/hr  filter  tuning  overestimates 
errors  grossly,  compared  with  the  0.4  or  2.0  nm/hr  filters.  Overestimating  of  the  true 
error  by  the  system  filter  will  be  discussed  in  Chapter  5. 
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4.8  Performance  Example  of  a  4.0  nm/hr  INS  during  Two  Periods  of  GPS 
Measurement  Loss  with  Radar  Altimeter  Aiding  during  the  Final  Approach. 

•  Case  XVI  (Figures  U  -  I  -U  -  6  in  Appendix  U) 

This  is  the  final  performance  evaluation  researched  in  this  thesis.  The  following  scenario 
take  place  using  the  Tanker  profile  (See  Section  4.1): 

•  T=  1000  seconds:  All  GPS  measurements  are  lost. 

•  T=1 120  seconds:  All  GPS  measurements  have  been  re-acquired.  A  120  second 
outage  was  chosen  to  "deprive"  the  INS  of  120  SV  measurements.  It  was  thought  this 
was  could  simulate  a  "popped”  GPS  receiver  circuit  breaker  on  the  aircraft,  allowing 
the  GPS  receiver  time  to  re-acquire  the  same  SV  set  as  used  prior  to  losing  power. 

•  T=3300  All  GPS  measurement  are  lost. 

Aircraft  altitude  is  7524  feet  (MSL) 

Aircraft  altitude  is  6699  feet  (AGL) 

•  T=3641  seconds:  Radar  altimeter  measurements  are  available. 

Aircraft  altitude  is  4028  feet  (MSL) 

Aircraft  altitude  is  3203  feet  (AGL) 

•  T=3750  seconds:  GPS  measurements  are  available.  A  450  second  GPS  outage  was 
chosen  to  allow  enough  time  to  degrade  the  INS  performance  severely. 

Aircraft  altitude  is  2723  feet  (MSL) 

Aircraft  altitude  is  1898  feet  (AGL) 

•  T=3  9 1 5  Aircraft  lands  safely 

Aircraft  altitude  is  825  feet  (MSL) 

Aircraft  altitude  is  0  feet  (AGL) 

Table  4-23  summarizes  Case  XVI  landing  system  performance  at  respective  decision 
heights.  Table  4-24  summarizes  the  Case  XVI  performance  in  terms  of  what  category  of 
precision  approach  is  met  based  on  Horizontal  and  Vertical  errors. 


Case  XVI  was  compared  with  Case  VI.  Results  are  shown  in  Table  4  -  25.  Looking  at 
Table  4  -  25,  it  is  evident  that  a  Category  I  landing  is  feasible,  given  the  GPS  outages 
listed  above,  as  long  as  an  independent  measurement  -  the  radar  altimeter,  is  functioning. 
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Position  Error 
State 

Decision 
Height  (feet) 

Case  XVI  True 
Error  (feet) 

Case  XVI 
Filter  Error 
(feet) 

Slatitude 

200 

9.29 

9.78 

Slatitude 

100 

B.86 

9.78 

Slatitude 

50 

9.47 

9.78 

Slongitude 

200 

9.75 

10.7 

Slongitude 

100 

9.55 

10.7 

Slongitude 

50 

9.71 

10.7 

Saltitude 

200 

3.04 

11.51 

Saltitude 

100 

3.69 

11.4 

Saltitude 

50 

3.98 

11.4 

Table  4-23  Case  XVI  la  Latitude,  Longitude  and  Altitude  Errors 


Case 

Number 

Horizontal 
(Based  on  True 
Estimated  Error) 

Horizontal 
(Based  on  Filter 
Predicted  Error) 

Vertical 
(Based  on  True 
Estimated  Error) 

Vertical 

(Based  on  Filter 
Predicted  Error) 

XVI 

Cat  I 

None 

Call 

None 

Table  4  -24.  Precision  Landing  Category  Predicted  To  Achieve 


The  plots  of  Appendix  U  also  show  that  the  4.0  nm/hr  filter  tuning  overestimates  errors 
grossly,  compared  with  the  0.4  or  2.0  nm/hr  filters.  Overestimating  of  the  true  error  by 
the  system  filter  will  be  discussed  in  Chapter  5. 


Error  State 

Case  VI  average  true 
error  (feet) 

Case  XVI  average  true 
error  (feet) 

%  Change  in  Error 

Latitude 

8.8 

9.2 

4.5% 

Longitude 

9.4 

9.7 

3.2% 

Altitude 

2.3 

3.6 

56.5% 

Clock  Bias 

6.0 

6.1 

1.7% 

Table  4  -  25.  Averaged  True  Error  %  Change,  Case  VI  vs.  Case  XVI 


4.9  Chapter  Summary 


This  ends  the  analysis  of  the  data  collected  from  Case  I  to  Case  XVI  scenarios.  Tables  4  - 
26  and  4-27  summarize  the  results  from  Case  I  -  XVI,  respectively.  Note  that  no 
Category  IH  precision  approach  was  deemed  possible  by  any  Case  number  in  this  thesis. 
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The  next  chapter  summarizes  the  results  of  this  thesis  and  presents  recommendations  for 
future  research. 


Precision 

Approach 

Category 

Case 

I 

Case 

II 

Case 

III 

Case 

IV 

Case 

V 

Case 

VI 

Case 

VII 

Case 

VIII 

I 

No 

Yes 

No 

Yes 

No 

Yes 

No 

No 

II 

i 

1 

i 

— 

— 

— 

i 

i 

i 

i 

l 

i 

— 

i 

i 

Table  4  -  26.  Summary  of  Cases  I  -  VIII:  Precision  Approach  Requirements  Met 


Case 

IX 

Case 

X 

Case 

XI 

Case 

XII 

Case 

XIII 

Case 

XIV 

Case 

XV 

Case 

XVI 

I 

No 

_ 

l 

l 

No 

No 

No 

Yes 

II 

i 

i 

i 

Yes 

Yes 

Yes 

— 

i 

i 

i 

i 

i 

i 

l 

l 

l 

Table  4  -  27.  Summary  of  Cases  IX  -  XVI:  Precision  Approach  Requirements  Met 
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V.  Conclusions  and  Recommendations 


Chapter  5  summarizes  the  research  and  presents  recommendations  and  conclusions  based 
on  results  presented  in  Chapter  4. 

5.1  Introduction 

In  the  case  of  aircraft  integration,  very  few  "integrated"  GPS/INS  studies  have 
been  accomplished  with  respect  to  implementing  precision  landing  approaches.  Instead, 
the  majority  of  research  efforts  explored  stand-alone  GPS  receiver  technology.  Instead  of 
stand-alone  GPS  techniques,  this  thesis  integrated  the  GPS  with  an  INS/Baro  altimeter 
system  and  a  radar  altimeter  using  an  extended  Kalman  filter  to  meet  FAA  Category  I  and 
II  precision  approach  accuracy  requirements. 

In  order  to  accomplish  thesis  goals,  this  project  concentrated  on  setting  up  reliable 
models  for  three  different  classes  of  INS  (0.4,  2.0  and  4.0  nm/hr  (CEP)  INS  systems), 

GPS  (4  channel  receiver,  5m  vertical  and  4m  horizontal  precision,  la),  Baro  Altimeter 
(50  -  150  ft,  la),  and  Radar  altimeter  (1%  of  altitude  +  "floor"  value,  la).  This  thesis 
also  developed  a  generic  precision  approach  flight  profile  (using  PROFGEN  [47])  that 
encompassed  a  majority  of  aircraft  types.  Lastly,  this  thesis  utilized  a  single  ground- 
based  SV  (pseudolite)  and  available  true  post-processed  ephemeris  data  [14,24],  instead 
of  prior  simulated  ephemeris  data  used  at  AFIT  [58,3,45,50,27].  Once  all  the  above 
elements  were  in  place,  the  Multimode  Simulation  for  Optimal  Filter  Evaluation 
(MSOFE)  [48]  was  utilized  to  perform  extended  Kalman  Filter  integration  analysis. 

5.2  Conclusions 

Use  of  an  existing  0.4  nm/hr,  2.0  nm/hr  and  lower  quality  4.0  nm/hr  INS,  when  properly 
integrated  with  a  4-channel  P-code  GPS  receiver  and  radar  altimeter,  can  meet  the  FAA 
Category  I  precision  approach.  Use  of  an  existing  0.4  nm/hr,  2.0  nm/hr  and  lower  quality 
4.0  nm/hr  INS,  when  properly  integrated  with  a  5-channel  P-code  GPS  receiver  (one 
channel  using  a  ground  based  pseudolite)  and  radar  altimeter,  can  meet  the  FAA  Category 
II  precision  approach.  These  conclusions  are  made  based  on  mainly  four  main  factors: 

1.  Error  models  used  in  this  simulation  are  realistic  to  the  respective  real-world  black 
box  output  errors. 
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2.  No  radar  altimeter  measurement  outages  occur  during  the  landing  approach. 

3.  When  radar  altimeter  measurements  are  available,  the  earth’s  surface  will  be  modeled 
as  flat  and  referenced  approximately  to  the  INS-indicated  altitude  (referenced  to 
WGS-84  ellipsoid). 

4.  When  use  of  the  single  pseudolite  information  is  used,  ionospheric,  tropospheric, 
pseudolite  position,  pseudolite  clock  and  multipath  errors  are  negligible. 

The  use  of  real  ephemeris  data  obtained  from  the  National  Geodetic  Office  was  a  first  for 
use  in  an  AFIT  thesis.  It  is  recommended  by  the  author  to  use  real  ephemeris  data,  rather 
than  use  existing  FORTRAN  functions  that  simulate  ephemeris  data  when  possible. 

The  tanker  profile  looks  to  be  a  realistic  flight  scenario  of  a  precision  ILS  approach  at 
Wright-Patterson  AFB,  OH.  PROFGEN  [47]  was  utilized  to  create  a  trajectory  that 
would  represent  a  generic  transport  aircraft  from  take-off  to  a  precision  approach  landing. 
This  flight  profile  was  named  “Tanker”  since  its  attributes  closely  resemble  those  of  a 
KC-135  tanker  (Boeing  707)  aircraft. 

The  two  additional  INS  models  (2.0  nm/hr  and  4.0  nm/hr)  were  created  based  on  the  NRS 
model  of  [45].  All  models  represent  a  set  of  "generic"  INS's.  All  performance 
evaluations  conducted  in  this  research  showed  the  0.4,  2.0  and  4.0  nm/hr  INS  performed 
equally  well  until  loss  of  GPS  measurements  occurred.  If  loss  of  GPS  measurements 
occurred,  the  lower  quality  INS's  navigation  performance  degraded  rapidly.  Lastly,  due  to 
the  high  "Rfiiter"  radar  altimeter  measurement  noise  value  and  pseudolite  measurement 
noise  value,  performance  of  the  filter  was  overly  conservative  (achieving  la  errors  of  5  - 
7  feet)  during  the  final  landing  approach.  Proper  tuning  of  the  Rfiiter  measurement  noise 
values  would  allow  the  filter's  estimate  of  true  errors  to  be  less  conservative  (and  more 
realistic).  Time  permitting,  additional  tuning  of  the  filter  model  would  have  been 
accomplished  in  this  thesis. 

5.3  Recommendations 

This  thesis  has  pushed  the  performance  capabilities  of  the  hardware  stated  in  Section  5.2 
to  meet  a  Category  I  and  II  specifications  for  landing.  In  order  to  meet  a  Category  HI 
approach,  more  precise  measurements  must  be  made  available.  Recommendations  are  as 
follows: 
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1.  Make  use  of  the  centimeter  accuracy  of  carrier-phase  GPS  signals,  now  readily 
available  in  commercial  receivers.  The  carrier-phase  information  can  be  used  for 
positional  information  as  well  as  heading  information  (using  multiple  antennas)  [64]. 
Decide  whether  one's  algorithm  will  also  handle  "integer  ambiguity"  techniques,  or 
assume  no  cycle  slips.  See  Section  1.3.2  or  [27,6]  for  more  information  regarding 
carrier-phase  integer  ambiguity  and  cycle  slips. 

2.  Use  Differential  GPS  measurements  at  all  times  for  North  America  (assume  Wide- 
area  differential  will  be  available).  The  models  for  filter  design  and  for  performance 
evaluation  must  make  the  distinction  of  representing  a  C/A  or  P-code  receiver  using 
differential  corrections. 

3.  Maintain  use  of  at  least  one  pseudolite  along  the  flight  path.  In  fact,  current  flight 
testing  is  showing  that  use  of  two  pseudolites  optimally  placed  along  the  approach 
path  is  recommended  [11]. 

4.  Use,  as  a  minimum,  six  (6)  independent  SV  measurements  (though,  preferably,  for 
any  given  time,  use  of  all  in-view  SVs  is  preferred)  rather  than  the  current  (archaic) 
military  standard  of  four  SVs  (to  minimize  geometric  dilution  of  precision  (GDOP 
errors).  (This  thesis  only  used  cases  where  4  SV  or  5  SV  were  used  at  one  time). 

5.  Obtain  and  use  published  geographic  data  of  Wright-Patterson  AFB,  and  vicinity,  so 
as  to  reference  the  radar  altimeter  outputs  to  WGS-84  ellipsoid  (to  match  the  INS 
positional  outputs). 

6.  Use  single  filter  (then  continue  analysis  with  small,  non-computer-burdening  multiple 
filters)  and  perform  residual  monitoring  for  use  as  a  fault  detection  and  isolation 
algorithm  for  GPS  Space  Segment  system  errors,  to  compensate  for  deliberate  and 
non-deliberate  jamming  and  spoofing  of  the  SV  signals.  Fault  detection  must  notify 
the  pilot  of  a  possible  degraded  navigation  solution  in  less  than  2  seconds,  while 
minimizing  false  alarms. 

7.  Explore  utilizing  dynamic  filter  tuning  procedures  along  the  approach  path  when 
using  pseudolites  (i.e.,  at  a  given  runway,  multipath  errors  may  be  excessive;  "R" 
tuning  of  the  respective  measurement  may  be  necessary). 
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8.  All  analysis  should  show  a  detailed  window  of  performance  based  on  realistic  SV 
geometries.  SV  geometries  used  in  the  analysis  should  be  a  function  of  the  predicted 
number  of  SVs  available  within  +  5  years  time  (based  on  SV  launch  schedules). 
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Appendix  A.  Error  State  Definitions  for  the  LSM  Truth  and  Filter  Models 


Tabular  listings  of  the  truth  and  filter  models  are  presented.  Tables  A  - 1  and  A  -  2  show 
the  39  INS  states  for  the  truth  model,  with  the  SNU  84-1  [litton]  and  NRS  [mosle]  state 
numbers  given  for  cross-reference.  Table  A.  3  list  the  GPS  states  respectively,  and  Table 
A  -  4  lists  the  states  in  the  reduced-ordered  LSM  filter  model. 

Note:  In  Table  A  -  1,  the  LSM  states  12  and  13  are  not  included;  these  are  found  in  Table 
A -3. 
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LSM 

State 

State 

Symbol 

Definition 

SNU  84-1 
State 

NRS 

State 

1 

KB 

X-component  of  vector  angle  from  true  to  computer  frame 

1 

1 

2 

■a 

Y-component  of  vector  angle  from  true  to  computer  frame 

2 

2 

3 

Z-component  of  vector  angle  from  true  to  computer  frame 

3 

3 

4 

X-  component  of  vector  angle  from  true  to  platform  frame 

4 

4 

5 

Y-  component  of  vector  angle  from  true  to  platform  frame 

5 

5 

6 

Z-  component  of  vector  angle  from  true  to  platform  frame 

6 

6 

7 

Hsu 

X-component  of  error  in  computed  velocity 

7 

7 

1  8 

fiS 

Y-component  of  error  in  computed  velocity 

8 

8 

9 

mvm 

Z-component  of  error  in  computed  velocity 

9 

9 

10 

8h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

10 

10 

11 

She 

Total  baro-altimeter  correclated  error 

23 

11 

14 

ShL 

Error  in  lagged  inertial 

11 

16 

15 

msm 

Error  in  vertical  channel  aiding  state 

12 

17 

16 

8S4 

Error  in  vertical  channel  aiding  state 

13 

18 

17 

Axc 

X-component  of  accelerometer  and  velocity  quantizer 
correlated  noise 

17 

19 

18 

Ayc 

Y-component  of  accelerometer  and  velocity  quantizer 
correlated  noise 

18 

20 

19 

Azc 

Z-component  of  accelerometer  and  velocity  quantizer 
correlated  noise 

19 

21 

20 

E9 

X-component  of  gravity  vector  errors 

20 

22 

21 

Y-component  of  gravity  vector  errors 

21 

23 

22 

Sgz 

Z-component  of  gravity  vector  errors 

22 

24 

Table  A  - 1.  39-State  INS  System  Model:  First  20  States 


Note:  LSM  state  12  and  state  13  are  located  in  Table  A  -  3 


State 

Number 


Definition 


X-component  of  gyro  drift  rate  reoeatabili 


Y-component  of  gyro  drift  rate  repeatabili 


_  X-component  of  gyro  scale  factor  error 


_  Y-component  of  gyro  scale  factor  error 


Z-component  of  gyro  scale  factor  error 


X-component  of  accelerometer  bias  repeatability 


Y-component  of  accelerometer  bias  repeatability 


Z-component  of  accelerometer  bias  repeatability 


Sax  X-component  of  accelerometer  and  velocity  quantizer 

scale  factor  error 


SAv  Y-component  of  accelerometer  and  velocity  quantizer 

scale  factor  error 


Saz  Z-component  of  accelerometer  and  velocity  quantizer 

scale  factor  error 


S()Ax  X-component  of  accelerometer  and  velocity  quantizer 
_ scale  factor  asymmetry 


SoAv  Y-component  of  accelerometer  and  velocity  quantizer 
_ scale  factor  asymmetry 


SoAz  Z-component  of  accelerometer  and  velocity  quantizer 
_ scale  factor  asymmetry 


_  X  accelerometer  misalignment  about  Z-axis 


_ Y  accelerometer  misalignment  about  Z-axis 


Z  accelerometer  misalignment  about  Y-axis 


X  accelerometer  misalignment  about  Y-axis 


SNU  84-1 
State 

NRS 

State 

30 

25 

31 

26 

32 

27 

33 

28 

34 

29 

35 

30 

48 

31 

49 

32 

50 

33 

51 

34 

52 

35 

53 

36 

54 

37 

55 

38 

56 

39 

66 

40 

67 

41 

68 

42 

69 

43 

Table  A  -  2.  39-state  INS  System  Model:  Second  19  States 


A  -  3 


State 

Number 


State  Symbol 


Defiuition 


NRS  State 


8Rclku 


8 


SRcloonl 


SRtroo  1 


8Rjon  1 


SRclksvl 


8xSvl 


8ySvl 


8zsvl 


8 


8Rtroo2 


8Rjon2 


5Rc 


8 


SySv2 


8zSv2 


8 


3  _ 


SB 


8Rtrop4 


n4 


8Rclksv4 


8 


5ysv4 


User  clock  bias 

14 

User  clock  drift 

15 

SV  1  code  loop  error 

68 

SV  1  tropospheric  error 

69 

SV  1  ionospheric  error 

70 

SV  1  clock  error 

71 

SV  1  x-component  of  position  error 

72 

SV  1  y-component  of  position  error 

73 

SV  1  z-component  of  position  error 

74 

SV  2  code  loop  error 

75 

SV  2  tropospheric  error 

76 

SV  2  ionospheric  error 

77 

SV  2  clock  error 

78 

SV  2  x-component  of  position  error 

79 

SV  2  y-component  of  position  error 

80 

SV  2  z-component  of  position  error 

81 

SV  3  code  loop  error 

82 

SV  3  tropospheric  error 

83 

SV  3  ionospheric  error 

84 

SV  3  clock  error 

85 

SV  3  x-component  of  position  error 

86 

SV  3  y-component  of  position  error 

87 

SV  3  z-component  of  position  error 

88 

SV  4  code  loop  error 

89 

SV  4  tropospheric  error 

90 

SV  4  ionospheric  error 

91 

SV  4  clock  error 

92 

SV  4  x-component  of  position  error 

93 

SV  4  y-component  of  position  error 

94 

SV  4  z-component  of  position  error 

95 

State 

Number 

State 

Symbol 

Definition 

NRS  State 

1 

80x 

X-component  of  vector  angle  from  true  to  computer  frame 

1 

2 

50v 

Y-component  of  vector  angle  from  true  to  computer  frame 

2 

3 

807 

Z-component  of  vector  angle  from  true  to  computer  frame 

3 

4 

4>x 

X-  component  of  vector  angle  from  true  to  platform  frame 

4 

5 

Y-  component  of  vector  angle  from  true  to  platform  frame 

5 

6 

Z-  component  of  vector  angle  from  true  to  platform  frame 

6  | 

7 

■■ 

X-component  of  error  in  computed  velocity 

7 

8 

■Kszm 

Y-component  of  error  in  computed  velocity 

8 

9 

UKimi 

Z-component  of  error  in  computed  velocity 

9 

10 

8h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

10 

11 

She 

Total  baro-altimeter  correclated  error 

11 

12 

User  clock  bias 

14 

13 

■E21ZI9B 

User  clock  drift 

15 

Table  A  -  4.  13-State  Reduced-Order  Filter  Model 


Appendix  B.  Dynamics  Matrices  and  Noise  Values 


B.  1  Definition  of  Dynamics  Matrices 

In  Chapter  3,  the  truth  and  filter  model  dynamics  are  defined  by  the  submatrices, 

F Filter'  FiNstl  .  FiNSn  and  fgps,  of  Equation  (3.3).  The  F rater  represents  the  filter 

dynamics  matrix,  which  is  also  a  submatrix  of  the  larger  truth  model  dynamics  matrix 
[mosle].  The  other  three  matrices  represent  the  additional  truth  model  non-zero  portions 
of  the  F  matrix  that  simulate  the  real  world  [mosle].  Tables  B.2,  B.3,  B.4  and  B.5  contain 
the  non-zero  elements  of  the  dynamics  submatrices  FFuter,  FtNS„ »  fins,2  and  fgps,  * 

respectively.  All  undeclared  variables  shown  in  the  following  tables  are  defined  in  the 
LN-93  technical  report,  along  with  their  units  [litton, mosle].  The  structure  of  the 
dynamics  matrices  in  this  chapter  correspond  to  the  truth  model  state  definitions  in 
Appendix  A  and  to  the  AFIT  thesis  (NRS  model)  by  [mosle].  The  notation  used  in 
Tables  B.2  -  B.4  in  this  Appendix  is  defined  in  Table  B.l. 
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Px.  py,  pz 

Components  of  angular  rate  of  navigation  frame  with 
respect  to  the  earth  (craft  rate),  coordinatized  in  Litton 
True  frame 

Components  of  earth  sideral  rate  vector  (earth  rate), 
coordinatized  in  Litton  ECEF,  with  respect  to  inertial 

space 

a 

Equatorial  radius  of  the  earth  (6378388  meters) 

_ go _ 

Equatorial  gravity  magnitude  (32.08744  ft/sec^) 

^itx,  Cfyty,  ^itz 

Components  of  angular  rate  of  navigation  frame  with 
respect  to  inertial  space  (spatial  rate),  coordinatized  in 
Litton  True  frame 

vx.vv.vz 

Components  of  vehicle  velocity  vector  with  respect  to 
earth-fixed  coordinates 

az 

Components  of  specific  force,  coordinatized  in  Litton 
True  frame 

Crx,  Cry 

Components  of  earth  spheroid  inverse  radii  of  curvature 

®ibx,  ®iby,  ®ibz 

Components  of  angular  rate  of  navigation  frame  with 
respect  to  inertial  space  (spatial  rate),  coordinatized  in 
Litton  Body  frame 

Q.i 

Elements  of  the  transformation  matrix  C™vdy 

Pbhc 

Barometer  inverse  correlation  time  (600  seconds) 

PVxc,  PV yc,  PVzc 

Gyro  inverse  correlation  time  constants  (5  minutes) 

P8gx,  P8gy,  PSgz 

Gravity  vector  error  inverse  correlation  time  constants 
(Velocity/correlation  distance) 

_2 

CTS  he 

Variance  of  barometer  correlated  noise 

2  2  2 
^Vxc’  ®Vyc  »  ®Vzc 

Variances  of  accelerometer  correlated  noise 

2  2  2 
c8gx’  ®8gy’  °8gz 

Variances  of  gravity  vector  correlated  noise 

2  2  2 

CTilto’  CT\by’®r\bz 

Power  spectral  density  value  of  gyro  drift  rate  white 

noise 

Power  spectral  density  value  of  accelerometer  white 

noise 

ki,k2,  k3,  k4 

Vertical  channel  gains  of  vertical  channel  error  model 
(see  figure  2  of  [litton]) 

Table  B  -  1 .  Notation  of  Variables  used  in  Tables  B  -  2  to  B  -  4 


B  -  2 


B  -  3 


Element 

Variable 

Element 

(9,16) 

-k2 

(9,17) 

(10,9) 

1 

(16,10) 

1 

1HW 

(17,18) 

-k3 

(18,10) 

(7,19) 

c„ 

(7,20) 

(7,22) 

1 

(8,19) 

(8,21) 

C33 

(8,23) 

(9,20) 

C32 

(9,21) 

(9,11) 

k2 

(10,11) 

(18,11) 

k4/600 

(18,18) 

(4,25) 

c„ 

(4,26) 

(4,28) 

CuC0ibx 

(4,29) 

(5,25) 

C21 

(5,26) 

(5,28) 

C2i(0ibx 

(5,29) 

(6,25) 

C31 

(6,26) 

(6,28) 

Qi^ibx 

(6,29) 

(7,31) 

Cm 

(7,32) 

(7,34) 

Cu  Af 

(7,35) 

(7,37) 

■2 

M 

■ 

(7,38) 

(7,40) 

msaa$- 

(7,41) 

(7,43) 

Cl  3axs 

(8,31) 

(8,33) 

C23 

(8,34) 

(8,36) 

c  23  Af 

(8,37) 

(8,39) 

C23 

Af 

(8,40) 

(8,42) 

mm 

(8,43) 

(9,32) 

c32 

(9,33) 

(9,35) 

c32a^ 

(9,36) 

(9,38) 

ES 

■ 

(9,39) 

(9,41) 

-C32  A‘ 

(9,42) 

Table  B-3.  Elements  of  the 


B 


Variable 

Element 

Variable 

-1 

(9,18) 

k2 

(10,18) 

k,-l 

-1 

(17,16) 

k3 

k4 

(18,16) 

-k4 

C|2 

(7,21) 

C13 

C2, 

(8,20) 

C22 

1 

(9,19) 

C3| 

c33 

(9,24) 

1 

k, 

(17,11) 

-k3 

IC4-I 

(9,43) 

C33  Af 

Cl2 

(4,27) 

Ci3 

(4,30) 

Cl3(0ibz 

c22 

(5,27) 

C23 

(5,30) 

■3ns 

C32 

mwm 

C33  I 

(6,30) 

C12 

(7,33) 

C,3 

C,2  A“ 

(7,36) 

C,3Af 

m 

m 

■ 

(7,39) 

m 

El 

■ 

-c12a? 

(7,42) 

M3SSM 

C21 

(8,32) 

c22 

C21A® 

(8,35) 

C21 

Af 

(8,38) 

C22 

ab 

1 

c21a^ 

(8,41) 

-C22Af 

c23Af 

(9,31) 

C31 

C33 

(9,34) 

Ca.Af 

C33  A^/ 

(9,37) 

12 

1 

■ 

C33 

Af 

(9,40) 

C3iA* 

C33  A? 

— 

— 

;  Dynamics  Submatrix  FINSn 


-4 


Element 

Variable 

(19,19) 

-PVxc  =  -3.33E-3  sec-1 

(20,20) 

-PVyc  =  -3.33E-3  sec-1 

(21,21) 

-PVzc  =  -3.33E-3  sec-1 

(22,22) 

-P8gx  =  -8.22E-6 *tJVx  +  V2  +  V2  ft~l 

(23,23) 

-P5gy  =  -8.22E-6*V^2  +  V?  +  ^2  fi~l 

(24,24) 

-p5gz  =  -8.22E-6*V^2  +  ^2  +  ^2 

Table  B-4.  Elements  of  the  Dynamics  Submatrix  Fins, 


Element 

Variable 

Element 

Variable 

""" 

Element 

Variable 

(42,42) 

-1  ft2  /sec 

(43,43) 

-1/500  ft2  /  sec 

(44,44) 

(49,49) 

(50,50) 

- 1/500  ft2  /  sec 

(51,51) 

-1/1500  ft2/ sec 

(56,56) 

- 1  ft 2  /  sec 

(57,57) 

EJEEEBM 

■ 

(58,58) 

(63,63) 

nmn 

(64,64) 

■ 

(65,65) 

Table  B-5.  Elements  of  the  Dynamics  Submatrix  FGPS 


B.2  Elements  of  the  Process  Noise  and  Measurement  Noise  Matrices 

This  section  defines  the  dynamic  noise  strengths  and  measurement  noise  variances 
for  the  truth  and  filter  models.  The  truth  model  non-zero  dynamics  noise  strengths  are 
defined  in  Tables  B.6  and  B.7.  These  noise  strengths  correspond  to  the  driving  noises 
w Filter’ w ins,  and  wgps,  in  Equation  (3.3).  Note  that  the  a2  terms  in  Table  B.6  are 

variable  names  as  defined  in  the  Litton  technical  report  [litton]  and  do  not  represent 
variance  terms  typically  associated  with  the  notation  a2  [mosle].  The  filter  dynamics 
driving  noise  terms  implemented  after  filter  tuning  for  each  respective  INS  integration 
(0.4  nm/hr,  2.0  nm/hr  and  4.0  nm/hr)  are  listed  in  Tables  B.8,  B.9  and  B.10.  Finally,  the 
measurement  noise  variances  used  in  the  truth  and  filter  models  are  presented  in  Table 
B.l  1. 
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Element 

Variable 

(4,4) 

al  =  190.4E-15 

*1  bx 

[ft3  /  sec5] 

(7,7) 

oj  =  102.9E-9 
[ft2  /  sec4] 

(11,11) 

2P8fora8/ic= 
33.34  [ft2  / sec] 

(21,21) 

2P VzcaVzc  = 
2.75E-1 1 
[ft2  /  sec5] 

(24,24) 

2P8gz<78gz  = 
3.10E-13 
[ft3  / sec5] 

Element 

Variable 

(5,5) 

0^  =  190.4E-15 
[ft2  / sec4] 

(8,8) 

ai  =  102.9E-9 

M  Ay 

[ft2  /  sec4] 

(19,19) 

2PvxcaV;tr  = 
2.75E-1 1 
[ft2  /sec5] 

(22,22) 

2p8gxa8gx  = 
3.10E-13 
[ft3  / sec5] 

Element 

Variable 

(6,6) 

Cf^  =  190.4E-15 
]/r2  /  sec4] 

(9,9) 

ol  =  102.9E-9 
rUz 

[ft2  /  sec4] 

(20,20) 

2PvycaVyc  = 
2.75E-1 1 
[ft2  /  sec5] 

(23,23) 

^Psgy^Sgy  = 
3.10E-13 
[ft3  /sec5] 

Table  B-6.  Elements  of  Truth  Model  Process  Noise  Submatrix  for  the  INS  Truth  Model 


Element 

Variable 

Element 

Variable 

Element 

Variable 

(42,42) 

I1!^3^2b13 

(43,43) 

(44,44) 

0.004  ft 2  /  sec 

(49,49) 

j 

(50,50) 

(51,51) 

0.004  ft2/ sec 

(56,56) 

0.5  ft2  /  sec 

(57,57) 

vmmmB 

(58,58) 

(63,63) 

(64,64) 

(65,65) 

0.004  ft2 /see 

Table  B-7.  Elements  of  Trut 

a  Model  Process  N 

oise  for  GPS  States 

Table  B-8.  Filter  Process  Noise  Q  Values  for  Case  1  (using  0.4  nm/hr  INS) 


Element 


(11,11) 


(13,13) 


Variable 


5E-15  ft2/sec3 


Table  B-9.  Filter  Process  Noise  Q  Values  for  Case  2  (using  2.0  nm/hr  INS) 


Table  B-10.  Filter  Process  Noise  Q  Values  for  Case  3  (using  4.0  nm/hr  INS) 
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Measurement 

Truth  Noise 

Filter  Noise 

Baro  Altimeter 

2500  ft2 

3500  ft2 

Satellite  Vehicles 

2  ft2 

2 

75  ft 

Radar  Altimeter 

(See  function,  Chapter  3) 

(See  function,  Chapter  3) 

Table  B-l  1.  Truth  and  Filter  Measurement  Noise  R  Values  for  Cases  1-3 
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Appendix  C  PROF_IN  Input  File 


The  PROF_EN  input  file  for  PROFGEN  representing  the  Tanker  (Boeing  707)  Flight 
Profile  is  found  in  this  section. 


PR0F_IN 


Group  1,  General  Information 


This  is  a  PROFGEN  input  file  and  these  are  general  instructions  for  under¬ 
standing  its  organization  so  it  can  be  used  as  a  template  for  building 
similar  files. 

PROF_IN  is  the  parameter  control  file  for  PROFGEN.  This  file  is  self- 
documenting  and  consists  of  text  and  data  entered  in  list  directed  form. 
PR0F_IN  is  divided  into  three  Groups  with  contents  as  follows: 

Group  1:  General  Information; 

Group  2:  Problem  Control  Parameters; 

Group  3 :  Segment  Control  Parameters . 

Each  Group  is  required  and  each  is  constructed  in  two  parts.  The  first 
part  is  a  “REMINDER"  (such  as  this)  that  is  placed  at  the  head  of  the  Group 
to  identify  that  Group  and  convey  information  about  it.  The  second  part  is 
a  series  of  entries  for  the  parameters  in  that  Group.  Each  parameter  is 
defined  by  a  leading  "NOTE"  that  conveys  this  set  of  five  characteristics: 

NAME  TYPE  DESCRIPTION  UNITS  DEFAULT  VALUE. 

The  parameter  value  is  entered  following  its  NOTE.  All  data  in  file 
PROF_IN  are  entered  in  list  directed  form.  You  must  include  both  the  Group 
REMINDERS  and  the  parameter  NOTES  (or  simple  substitutes  such  as  ' 1 )  in 
your  own  PROF_IN  file  to  place-hold  this  material.  The  PROFGEN  problem 
title  and  the  compiler  choice  are  the  only  parameters  in  Group  1. 


' PTITLE  CH*80  PROFGEN  run  title  -  *  *  ' 

'  Tanker 3  Profile  (  all  outputs  set  to  support  INS/GPS  problem  in  MSOFE  ) ' 

' FOR_77  CH*8  Compiler.  Controls  calls  to  date,  time  and  CPU  timer 

routines.  Choices:  DEC,  LAHEY,  UNIX.  -  1  'UNIX' 


Group  2,  Problem  Control  Parameters 


Note  that  each  entry  in  this  Group  is  echoed  in  Table  1  of  PROF_OUT,  converted 
to  internal  computational  units  of  radians/feet/seconds,  placed  in  Namelist 
PRDATA,  and  echoed  again  in  converted  form  as  Namelist  output  in  Table  4. 

See  Table  3  of  PROF_OUT  for  help  in  specifying  IPRSET  and  IRTSET. 


NAME  TYPE  DESCRIPTION  UNITS  DEFAULT 


' Problem  specs 


NSEG  IN  number  of  flight  segments  -  0 

' TSTART  DP  initial  time  of  trajectory  sec  0.D0 

' LLMECH  IN  local  level  mechanization  index: 

1  alpha  wander;  2  constant  alpha; 

3  unipolar;  4  free  azimuth  -  2 

'RTOL  DP  relative  tolerance  for  integration  -  0.D0 

'ATOL  DP  absolute  tolerance  for  integration  -  0.D0 


27 

0 


1 

0.1D-11 

0.1D-11 


'Earth  parameters 


ESQ  DP  earth  eccentricity  squared  -  6 . 69437999013D-3 '  / 

'REQ  DP  earth  semimajor  axis  m  6.378137D6  '  / 


Figure  C  - 1.  PROFJN  Input  File  for  PROFGEN 
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'WEI 

DP 

earth  sidereal  rate 

r/s 

7 . 292115D-5 

'  / 

'GM 

DP 

earth  gravitational  constant 

m3  /s2 

3 . 986005D14 

'  / 

'GEE 

DP 

converts  PACC/TACC  to  accel  units 

m/s2 

9.8 

'  32. DO 

'  RHOSL 

DP 

air  density  at  sea  level 

kg/  m3 

1.225 

'  0 . 002377D0 

'  ZETA 

DP 

exponent  in  air  density  function, 
rho  =  RHOSL  *  exp(-ZETA  *  alt) 

1/m 

1.1385D-4 

'  3 . 47015D-20 

'Vehicle  specs 


ROLRAT 

DP 

maximum  roll  rate 

d/s 

1 .  DO 

'  45. DO 

'  ROLTC 

DP 

roll-axis  time  constant 

sec 

1 .  DO 

'  1 .  DO 

'  VMASS 

DP 

vehicle  mass 

kg 

1 .  DO 

'  15. DO 

'  RAREA 

DP 

reference  area  for  drag 

m2 

1 .  DO 

'  1 .  DO 

'  CDRAG 

DP 

coefficient  of  drag 

~ 

1 .  DO 

'  O.ODO 

' State 

mitial 

conditions 

VETO 

DP 

initial  earth  velocity  magnitude 

m/s 

0 .  DO 

'  / 

'  ROLLO 

DP 

initial  aircraft  roll  angle 

deg 

0 .  DO 

'  / 

'  PITCHO 

DP 

initial  aircraft  pitch  angle 

deg 

0 .  DO 

'  / 

•HEADO 

DP 

initial  ground  path  heading 

deg 

0 .  DO 

'  232. DO 

'ALPHAO 

DP 

initial  alpha  angle 

deg 

0  .DO 

'  45. DO 

'GLATO 

DP 

initial  geographic  latitude 

deg 

0 .  DO 

'  39.833D0 

'  TLONO 

DP 

initial  terrestrial  longitude 

deg 

0 .  DO 

'  -84.033D0 

'  CLONO 

DP 

initial  celestial  longitude 

deg 

0 .  DO 

’  / 

'ALTO 

DP 

initial  altitude 

m 

0 .  DO 

'  825. DO 

'I/O  controls 

JPRNT 

IN 

print  (PROF_OUT  file)  control 

- 

1 

'  1 

'  JRITE 

IN 

write  (FLIGHT  file)  control 

- 

0 

'  1 

'  JPLOT 

IN 

plot  (META  file)  control 

- 

0 

'  / 

'  SIUNIT 

LG 

I/O  in  SI  units  (else  English) 

- 

.T. 

'  .F. 

'XVU 

RL 

x  viewpoint  for  3-D  plot 

inches 

-7.5 

■  / 

'  YVU 

RL 

y  viewpoint  for  3-D  plot 

inches 

-7.5 

■  / 

'ZVU 

RL 

z  viewpoint  for  3-D  plot 

inches 

18. 

■  / 

'  DELR 

DP 

spread  in  roll  angle  outputs 

deg 

0 .  DO 

•  / 

'DELP 

DP 

spread  in  pitch  angle  outputs 

deg 

0 .  DO 

■  / 

'DELY 

DP 

spread  in  yaw  angle  outputs 

deg 

0 .  DO 

•  / 

■PSEP 

DP 

min  separation  in  printed  outputs 

sec 

0  .  DO 

■  / 

'  WSEP 

DP 

min  separation  in  written  outputs 

sec 

0 .  DO 

■  / 

' IPRSET 

INO 

indices  of  printed  variables 

nd 

52*-l 

1,2,3,40,43,5,7,8,9,14,16,17,18  / 

'  IRTSET  IN ( )  indices  of  written  variables 

1,2,3,40,43,5,7,8,9,14,16,17,18 


nd 


52*-l 


Group  3 ,  Segment  Control  Parameters 


Note  that  each  entry  in  this  Group  is  echoed  in  Table  2  of  PROF_OUT,  converted 
to  internal  computational  units  of  radians/feet/seconds,  placed  in  Namelist 
SGDATA,  and  echoed  again  in  converted  form  as  Namelist  output  in  Table  5. 


NAME  TYPE  DESCRIPTION  UNITS  DEFAULT 


' SEGLNT  DP ( )  time  interval  of  each  segment  sec  0.D0' 

3. DO  30. DO  160. DO  329. DO  5. DO  25. DO 


371. DO 

35. DO 

400. DO 

35. DO 

340. DO 

35. DO 

600. DO 

35. DO 

120. DO 

35.  DO 

372. DO 

70.  DO 

40. DO 

35. DO 

360. DO 

45.  DO 

104. DO 

35.  DO 

39. DO 

25.  DO 

242. DO  / 

Figure  C  - 1  (Continued).  PROF_IN  Input  File  for  PROFGEN 
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KTURN 


4 


IN()  maneuver  index:  1  vertical;  2  horizontal; 

3  jink;  4  straight;  5  roll;  6  free  fall 

4  4  14  12 
4  2  4  2  4 
2  3  2  1  1 
4  2  111 
2  12  12 
4  / 

' NPATH  INO  path  index:  1  great  circle;  2  rhumb  line  -  2 

/ 

’ TACC  DP ( )  maximum  centrifugal  acceleration  during  turn  gees  0.D0' 

0 .DO  0 . DO  0.21D0  0.D0  0.9D0  0.9D0 

0 . DO  0.9D0  0 . DO  0.9D0  0.D0 

0.9D0  0 . DO  0.9D0  0.311D0  0.9D0 

0 . DO  0.9D0  0.63D0  0.9DO  0.9DO 

0.9D0  0.9D0  0.9D0  0.24D0  0.9D0 

0 .  DO  / 

'  PACC  DP ( )  signed  value  of  accel  along  velocity  vector  gees  0 . DO ' 

0 . DO  0.262D0  0.038D0  0.01D0  -O.0875D0  0.01D0 

0.01D0  -O.OODO  0.008D0  -0.0D0  0.D0 

-0.0D0  0 . DO  -O.OODO  -O.ODO  0.D0 

0.003D0  -0.125D0  -0.085D0  0.D0  0.015D0 

-O.OSDO  0 . DO  -0.05D0  -0.1882D0  -O.ODO 

-O.ODO  / 

' DELHED  DP()  for  horizontal  turns,  desired  change  in  heading 

angle;  for  jinking  maneuvers,  maximum  variation 
of  heading  angle  deg  0.D0' 

0 . DO  0 . DO  0 . DO  0 . DO  0.D0  45. DO 

0 .  DO  45.  DO  0 .  DO  45.  DO  0.D0 

45. DO  2. DO  45. DO  0.D0  0.D0 

0 . DO  90. DO  0 . DO  0.D0  0.D0 

-45. DO  0 . DO  45. DO  O.DO  45. DO 

0 .  DO  / 

'DELPIT  DP ( )  for  vertical  turns,  desired  change  in  pitch; 

for  j inking  maneuvers ,  period  of  maneuver  deg  or  sec  O.DO' 

O.DO  O.DO  5. DO  O.DO  -5. DO  O.DO 

O.DO  O.DO  O.DO  O.DO  O.DO 

O.DO  200. DO  O.DO  -5. DO  5. DO 

O.DO  O.DO  -5.0D0  5.0D0  -1.75D0 

O.DO  1.75D0  O.DO  -3. DO  O.DO 

O.DO  / 

' DELROL  DP ( )  desired  change  in  roll  angle  for  roll  maneuver  deg  O.DO' 

O.DO  O.DO  O.DO  O.DO  O.DO  O.DO 

O.DO  1 . DO  O.DO  1 . DO  O.DO 

O.DO  O.DO  O.DO  O.DO  O.DO 

O.DO  O.DO  O.DO  O.DO  O.DO 

1 . DO  O.DO  1 . DO  O.DO  O.DO 

O.DO  / 

' DTPRNT  DP ( )  time  interval  for  formatted  printing  on  PROF_OUT  sec  1.D8’ 

27*1. DO/ 

'DTKITE  DP ( )  time  interval  for  unformatted  writes  on  FLIGHT  sec  1.D8' 

27*1. DO  / 

'DTPLOT  DP ( )  time  interval  for  plotted  output  on  file  META  sec  1.D8’ 

27*4. DO  / 

' ICREDO  LG ( )  reset  kinematic  state  to  its  T START  value  -  .F.  ' 

/ 


Figure  C  - 1  (Continued).  PROF_IN  Input  File  for  PROFGEN 
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Appendix  D 


The  Tanker  Right  Profile  Plots  are  found  in  this  section. 
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xIO4 


Figure  D  - 1.  Lattitude,  Longitiude  and  Altitude  of  the  Tanker  Profile 
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(knots) 


0  500  1000  1500  2000  2500  3000  3500  4000 

Time  (sec) 


Figure  D  -  2.  3-D  Position,  Wander  Azimuth  and  Velocity  of  the  Tanker  Profile 
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Y-Velocity  (knots)  X-Velocity  (knots) 
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Figure  D  -  4.  Roll,  Pitch  and  Heading  Changes  of  the  Tanker  Profile 
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Appendix  E 


FORTRAN  source  code  ADDSV.for,  Sample  True  ephemeris  data  and  Almanac  file 
"051994A.AL3"  are  contained  in  this  section.  In  order  to  interpret  Figure  E  -2,  use  the 
following  template  shown  in  Table  E  - 1.  See  SEM3.6  User  Manual  for  the  template  of 
Figure  E  -  3. 


* 

1994  5 

21 

4  15 

0.00000000 

p 

6 

-3808.264821 

-26030.699349  -2199.205556 

179.659571 

* 

YEAR  MONTH 

DAY 

SV_HOURS  SV_MINUTES 

SV_SECONDS 

POSmONJNFO 

sv# 

ECEF_X(  in  km)  ECEF_Y  (in  km)  ECEF_Z  (in  km) 

SV  Clock  Offset 

(in  microsec) 

1994 

- 

YEAR 

5 

= 

MONTH 

21 

= 

DAY 

4 

= 

SV_HOURS 

15 

= 

SV_MINUTES 

0.00000000 

= 

SV_SECONDS 

P 

= 

POSITION JNFO 

6 

= 

SV# 

-3808.264821 

= 

ECEF_X(  in  km) 

-26030.699349 

= 

ECEF_Y  (in  km) 

-2199.205556 

= 

ECEF_Z  (in  km) 

179.659571 

= 

SV  Clock  Offset  (in  microsec) 

Table  E  -  1.  Template  for  Understanding  Figure  E  -  2. 
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nnn  no  nonnnnn  nononn 


*DECK  ADDSV 

PROGRAM  ADDSV 
C 

C  This  program  will  read  in  TANKER3  FLIGHT  profile,  and  merge 
C  the  information  with  real  SV  ephemeris  data  ovtained  from  the 
C  Coast  Guard  BBS  (National  Geodetic  Survey  Ephemeris) 

C 

C  The  output  file  created  will  be  called  "FLIGHT_TANKER3_21my94p4'' 
C  This  file  will  need  to  be  re-named  "FLIGHT"  so  MSOFE  will 
C  read  it  in  properly.  Put  the  "FLIGHT”  file  in  the  "runs" 

C  directory  with  your  MS0FE_1N  file. 

C 

C  Adopted  from  J  Solomon  add.for  code:  9  Sep  94  R.A.  Gray 

C 

C 

INTEGER  I,J,K 
REAL  EPHPOS(10000,12) 

CHARACTER  PDATE*10  ,  PTIME*10 
CHARACTER  PTITLE  *80,  TEMP  *80 
INTEGER*4  NRT,  NYT 
INTEGER*4  IDFLT(17),  IDCHEK(29) 

C 

DOUBLE  PRECISION  TKNOT(IOOOO),  UKNOT(  10000,29) 

C 

DATA  IDCHEK  /  1,  2,  3,  40,  40,  40,  43,  43,  43,  5,  7,  8,  9,  14, 

&  16,17,18,60,61,62,60,61,62,60, 

&  61,62,60,61,62/ 

DATA  NYT  729/ 


#OPEN  FILES 


OPEN  (  UNIT  =  20, 

&  FILE  =  FLIGHT, 

&  FORM  =  UNFORMATTED', 

&  STATUS  =  'OLD') 

OPEN  (UNIT  =21, 

&  FILE  =  '21my94p4', 

&  FORM  =  FORMATTED’, 

&  STATUS  =  'OLD') 

OPEN  ( UNIT  =  22, 

&  FILE  =  'FLIGHT  jrANKER3_21my94p4’, 

&  FORM  =  'UNFORMATTED',  &  STATUS  =  UNKNOWN’) 


#READ  FLIGHT  HEADER 


READ  (20)  PDATE,  PTIME,  PTITLE 
READ  (20)  NRT,  (IDFLT(I),I=1,NRT) 


#READ  FLIGHT  DATA 


DO  100  K=  1,3920 

READ  (20)  TKNOT(K),  (UKNOT(K,  J),  J  =  1,NRT) 

Figure  E  -  1.  FORTRAN  source  code  "ADDSV.for" 


E-2 


100  CONTINUE 
C 

C - 

C  #READ  GPS  DATA 

C - 

C 

DO  200 1  =  1,3920 

READ(21,*)  TEMP 

READ(21,201)EPHPOS(I,1),EPHPOS(I,2),EPHPOS(I,3) 

READ(21,201)  EPHPOS(I,4),EPHPOS(I,5),EPHPOS(I,6) 
READ(21,201)  EPHPOS(I,7),EPHPOS(I,8),EPHPOS(I,9) 
READ(21,201)  EPHPOS(I,10),EPHPOS(I,1 1), 

&  EPHPOS(I,12) 

C 

C  CONVERT  FROM  KM  TO  FEET, . . 

C 

EPHPOS(I,l)=  EPHPOS(I,l)  *  3.2801  *  1000 
EPHPOS(I,2)=  EPHPOS(l,2)  *  3.2801  *  1000 
EPHPOS(l,3)=  EPHPOS(l,3)  *  3.2801  *  1000 
EPHPOS(I,4)=  EPHPOS(l,4)  *  3.2801  *  1000 
EPHPOS(I,5)=  EPHPOS(I,5)  *  3.2801  *  1000 
EPHPOS(I,6)=  EPHPOS(I,6)  *  3.2801  *  1000 
EPHPOS(I,7)=  EPHPOS(I,7)  *  3.2801  *  1000 
EPHPOS(I,8)=  EPHPOS(I,8)  *  3.2801  *  1000 
EPHPOS(I,9)=  EPHPOS(I,9)  *  3.2801  *  1000 
EPHPOS(I,10)=  EPHPOS(I.IO)  *  3.2801  *  1000 
EPHPOS(1, 11)=  EPHPOS(1, 11)  *  3.2801  *  1000 
EPHPOS(I,12)=  EPHPOS(I,12)  *  3.2801  *  1000 

200  CONTINUE 
C 

201  FORMAT(5X,F13.6, 1  X,F1 3.6, 1  X,F1 3.6) 

C 

C - 

C  #ADD  GPS  DATA  TO  FLIGHT  DATA 

C - 

C 

DO  300 1  =  1,3920 
DO  300  J  =  1,12 

UKNOT(I,J+17)  =  DBLE(  EPHPOS(I,J) ) 

300  CONTINUE 

C 

C 

C - 

C  #OUTPUT  GPS  /  FLIGHT  DATA 
C- - 

c 

WRITE  (22)  PDATE,  PTIME,  PTITLE 
WRITE  (22)  NYT,  (IDCHEK(I),I=1,NYT) 

C 

DO  400  K=  1,3920 

WRITE(22)  TKNOT(K),  (UKNOT(K,  J),  J  =  1.NYT) 

400  CONTINUE 
C 

CLOSE(20),  CLOSE(21),  CLOSE(22) 

STOP  '  ADD  DONE.’ 

C 

END 


Figure  E  -  1  (Continued).  FORTRAN  source  code  "ADDSV.for" 


E-3 


*  1994  5  21  4  15  0.00000000 

P  6  -3808.264821  -26030.699349  -2199.205556  179.659571 

P  16  22228.030253  -2914.054829  14269.400710  -72.472353 
P  17  -1978.539062-16483.210409  20480.389383  -51.084448 
P  28 -14135.480147  5206.630899  21842.559082  14.685241 

*  1994  5  21  4  15  1.0000000 

P  6  -3807.902060  -26030.482798  -2202.379803  179.659588 

P  16  22229.666889  -2913.256477  14267.009251  -72.472371 

P  17  -1976.416118-16484.803199  20479.313405  -51.084449 
P  28 -14136.128577  5203.928088  21842.759458  14.685243 

*  1994  5  21  4  15  2.00000000 

P  6  -3807.539205  -26030.265875  -2205.554003  179.659604 

P  16  22231.303287  -2912.458317  14264.617490  -72.472389 
P  17  -1974.293373-16486.396028  20478.236982  -51.084450 
P  28 -14136.777174  5201.225288  21842.959367  14.685245 

*  1994  5  21  4  15  3.00000000 

P  6  -3807.176256-26030.048580  -2208.728154  179.659621 

P  16  22232.939448  -2911.660349  14262.225425  -72.472406 
p  17  -1972.170829-16487.988896  20477.160114  -51.084451 
P  28  -14137.425939  5198.522499  21843.158810  14.685247 

*  1994  5  21  4  15  4.00000000 

P  6  -3806.813214-26029.830911  -2211.902258  179.659638 

P  16  22234.575372  -2910.862573  14259.833058  -72.472424 
P  17  -1970.048484-16489.581803  20476.082800  -51.084452 
P  28 -14138.074872  5195.819722  21843.357787  14.685249 

*  1994  5  21  4  15  5.00000000 

P  6  -3806.450077  -26029.612871  -2215.076314  179.659654 

P  16  22236.211057  -2910.064990  14257.440387  -72.472442 
P  17  -1967.926339-16491.174748  20475.005041  -51.084453 

P  28 -14138.723973  5193.116956  21843.556298  14.685251 

*  1994  5  21  4  15  6.00000000 

P  6  -3806.086847-26029.394457  -2218.250322  179.659671 

P  16  22237.846505  -2909.267599  14255.047414  -72.472460 
P  17  -1965.804394-16492.767731  20473.926836  -51.084454 
P  28  -14139.373241  5190.414202  21843.754343  14.685253 

Figure  E  -  2.  Sample  from  National  Geodetic  OfficeTrue  Ephemeris  Data 


E-4 


25  051994A. AL3 

750  32768 


4. 03787490677082 E- 0003  -2 . 473 82546733 128E-0009 
-2.18499073561888E-0001  -3 . 81680029160470E-0001 
-2.47955322270000E-0005  0 . 00000000000000E+0000 


1 

32 

7 

3 .618717193 60000E-0003 
5.15368505900000E+0003 
9.34495603829844E-0001 

0 

1 

2 

13 

7 

1.28645896910000E-0002 

5.15360546900000E+0003 

-2.11397536351559E-0001 

0 

1 

4 

34 
7 

3 . 16715240480000E-0003 
5 . 15349951200000E+0003 
5. 56436181933 807E-0001 

0 

1 

5 

35 
7 

2.04277038570000E-0003 
5 . 15354834000000E+0003 
-8.79619968280586E-0001 
0 
1 

6 

36 
7 

6. 088733673 10000E-0003 
5 . 15367529300000E+0003 
8 . 0 6 14 2 54 073 67 4 8E- 0001 

0 

1 

7 

37 
7 

6 . 2 0 93 7 3 47 4 100 00E- 00 03 
5 . 15368798800000E+0003 
-5. 47001490125 07 6E- 0001 
0 
1 

9 

39 

7 

2. 6097297668000 0E-0003 
5 . 15371142600000E+0003 
9 . 39667574686501E-0001 

0 

1 


3.63542500182467E-0003 

4.34110874460587E-0001 

-9.05990600590000E-0005 


6.62805331970454E-0003 

-8.88224868033215E-0001 

2.76565551760000E-0005 


4 . 40408 688739 476E-0003 
4 . 38525794771209E-0001 
4 . 19 6 1669922 000OE-O005 


5. 2 8 14 62 6095 93 43 E- 0003 
7 . 8 44 84 77 03 92 502E- 0001 
1 . 81198120120000E-0004 


6. 33 81307947032 9E-0003 
7.72793384134664E-0001 
6. 9618225098000 0E-0004 


3 . 08419082134277E-0003 
1.10435722418997E-0001 
-1.81198120120000E-0005 


-2. 64844865770802 E- 0009 
-8.51847509087840E-0001 
-3 . 63797880710000E-0012 


-2.54658524778148E-0009 
-3. 98996607241964 E- 0001 
O.OOOOOOOOOOOOOOE+OOOO 


-2.64844865770802E-0009 
-7 . 478 85202727 681E-0001 
3 . 637 97 8 8 07 100 OOE- 00 12 


-2.48837742344679E-0009 

-9.51250588590547E-0001 

1.81898940350000E-0011 


-2 . 4774 634563 5220E-0009 
-8. 64231552535642 E-0001 
0 . OOOOOOOOOOOOOOE+OOOO 


-2.61934474550884E-0009 

-2.02198618779358E-0001 

O.OOOOOOOOOOOOOOE+OOOO 


12 

10 

3 

1. 450777053 80000E-0002  4 . 65393139879435E-0002  -2 . 19006317174084E-0009 


Figure  E  -  3.  SEM3.6  Almanac  Data  File  (051994.al3)  Used  For:  21  May  94. 


E-5 


5 . 15351025400000E+0003  2 . 20944414025568E-0001  -5 . 63353322948398E-0002 

6 . 88273 66592497 IE-0001  3 . 2424926758000QE-0005  0 . OOOOOOOOOOOOOOE+OOOO 


3.12614440920000E-0003 
5.15373291000000E+0003 
-6. 4970753 5115492E-0001 
0 
1 


5 .97 19188140793 4E-0003  -2 . 55022323 680240E-0009 
-5.47894613209864E-0001  9 . 780193 89301773E-0001 
4.76837158200000E-0006  0 . OOOOOOOOOOOOOOE+OOOO 


6 . 81734085080000E-0003 
5 . 15363085900000E+0003 
1 . 15476012945291E-0001 


1.10052823333991E-0003  -2 . 52111932463505E-0009 
1 .77  898912451352E-0001  5 . 72  649608752327E-0001 
!  .48770141600000E-0005  3 . 63797880710000E-0012 


8 . 64505767820000 E- 0004 


5 . 15353662100000E+0003  -5 . 44454702377718E-0001 
1.043 87400349130E-0001  -7 . 24792480470000E-0005 


4.97436971769650E-0003  -2 . 56841289921599E-0009 
5.44454702377718E-0001  -5 . 09196621865252E-0001 


0. OOOOOOOOOOOOOOE+OOOO 


7.39622116090000E-0003 
5 . 15364502000000E+0003 
8 . 3 05424 80168 52 OE-OOOl 


8 . 47055372544914E-0003  -2 . 51748133561413E-0009 
-8. 669319056683 87 E-0001  6 . 00263000159655E-0001 
-5.14984130860000E-0005  0 . OOOOOOOOOOOOOOE+OOOO 


5. 53703308110000 E-0003 
5.15374414100000E+0003 
-5.89919703533139E-0001 
0 
1 


3. 05195604 6725 14E-0 004  -2 . 51384334659322E-0009 
-2.32049219860062E-0001  4 . 18790684681791E-0001 
-7 . 62939453120000E-0006  0 . OOOOOOOOOOOOOOE+OOOO 


1 . 25408172610000E-0004 
5 . 15359960900000E+0003 
3.50051511019709E-0001 


-2. 17437168560882 E-0003  -2 . 66300061379169E-0009 
1.02508064366804E-0001  -8 . 89521010320165E-0001 
2.86102294920000E-0005  3 . 63797880710000E-0012 


4 . 7 693 252563 OOOOE-O003 
5 . 15353125000000E+0003 
-2.63615848922248E-0001 
0 


5. 535128239003 81E-0003  -2 . 63389670162434E-0009 
4.36904307926814E-0001  4 . 58547696903478E-0001 
5 . 34057617190000E-0005  0 . OOOOOOOOOOOOOOE+OOOO 


Figure  E  -  3  (Cont'd).  SEM3.6  Almanac  Data  File  (051994.al3)  Used  For:  21  May  94. 


E-6 


21 

21 

7 

1 . 103 9733 8870000E-0002 
5. 15362597700000E+0003 
9. 14788257969335 E- 0003 

0 

1 

22 

22 

7 

7.18355178830000E-0003 

5.15354003900000E+0003 

-2.24944476253860E-0001 

0 

1 

23 

23 
7 

8.32128524780000E-0003 

5.15359179700000E+0003 

-1.65312885266793E-0001 

0 

1 

24 

24 
7 

5.48076629640000E-0003 
5. 15362744100000E+0003 
6. 7847416613 6823E-0001 

0 

1 

25 

25 
7 

5.81312179570000E-0003 
5.15355957000000E+0003 
-7 . 64856704540744E-0001 
0 
1 

26 
26 
7 

8 . 3 184242249 OOOOE-0003 
5.15363476600000E+0003 
-4. 9720322733383 5E-0001 
0 
1 

27 

27 

7 

1 . 09333 992000000E-0002 
5 . 15369824200000E+0003 
5.49562924837777E-0001 

0 

1 


3.97873687860613E-0003 
-5. 57345393927682 E-0001 
-2.57492065430000E-0005 


3. 68310352 42 6479E-0003 
4.38924791650606E-0001 
1. 2397766113 0000E-0004 


5.13651081681166E-0003 
-5 . 46834947341484E-0001 
5.72204589840000E-0006 


9.78280379726186E-0003 
-8. 913991604 13 578 E-0001 
6 . 27 5177002 0000OE-O004 


4 . 9553 0209 660800E-0003 
-2.2  542 858487 1805E-0001 
-6.48498535160000E-0005 


2 . 11908510906553E-0003 
1. 0645735552 6676E-0001 
2.95639038090000E-0005 


-2.57568887725782E-0009 
8. 95620672 682456E-0001 
O.OOOOOOOOOOOOOOE+OOOO 


-2.65936262477077E-0009 
-7.99311402531010E-0002 
3 . 63797880710000E-0012 


-2 . 55749921484423E-0009 
-7.51196023020909E-0001 
O.OOOOOOOOOOOOOOE+OOOO 


-2.51020535757230E-0009 
-6. 87911362216682 E-0001 
3 . 27418092640000E-0011 


-2. 4556358 04 94953 E- 0009 
-3.30814701228596E-0001 
-3 . 6379788 07 10000E- 0012 


-2.61934474550884E-0009 

7.64407505673693E-0001 

O.OOOOOOOOOOOOOOE+OOOO 


8.94565027331417E-0004  -2 . 63753469064526E-0009 
1 . 02575664265170E-0001  9 . 12953036069076E-0001 

-7. 62939453 120000E-0006  O.OOOOOOOOOOOOOOE+OOOO 


28 

28 

7 


Figure  E  -  3  (Cont'd).  SEM3.6  Almanac  Data  File  (051994.al3)  Used  For:  21  May  94. 


E-7 


5 .21230697630000E-0003 
5 . 153 62744 100000E+0003 
3 . 86364096897161E-0001 

0 

1 

29 

29 

7 

5.16462326050000E-0003 
5.15363330100000E+0003 
5. 61828388466135 E- 0001 

0 

1 

31 

31 

7 

5.00249862670000E-0003 
5.15360058600000E+0003 
9. 34471242810989 E- 0001 

0 

1 


8. 884444196682 18E-0003 
7.75170821127847E-0001 
1.52587890620000E-0005 


3. 66022239169097 E- 0003 
-2.32212423429801E-0001 
1. 33514404300000 E- 0005 


6. 14358422218893 E-0003 
7.73816318762073E-0001 
1.52587890620000E-0005 


2 .45199781592861E-0009 
9.17673683276588E-0001 
3 . 63797880710000E-0012 


2. 473 8254673 3128E-0009 
5.82052691976437E-0001 
O.OOOOOOOOOOOOOOE+OOOO 


2. 481101445373 12 E- 0009 
2. 03706971195042 E- 0001 
3. 6379788 07 10000E- 0012 


Figure  E  -  3  (Cont'd).  SEM3.6  Almanac  Data  File  (051 994.al3)  Used  For:  21  May  94. 


Appendix  F 


Plots  obtained  from  SEM  3.6  (GDOP,  PDOP,  VDOP,  TDOP,  SV  Bearing/Elevation,  SV 
Rise/Set,  Number  of  Visible  SVs  and  Elevation/Time)  are  shown  in  this  section  (Figures 
F  -  1  to  F  -  5).  The  best  four  SV  selected  for  use  in  this  thesis  are  shown  in  Section  4.2. 
Notes:  In  Figure  F  -  3,  SV  Bearing/Elevation  Plot,  the  trajectory  starts  at  the  respective 
SV  number.  In  Figure  F  -  4  ,  SV  Rise/Set,  the  solid  line  shows  the  best  4  SV  for  a  given 
instant  in  time.  In  Figure  F  -  5,  Elevation/Time,  the  SV  number  lies  directly  under  the 
start  of  its  elevation  line.  For  more  information,  Use  reference  [18]. 


•  Used  in  Nav  Join 

*  Not  Used  in  Nav  Join 


Figure  F  -  3.  TDOP  and  SV  Bearing/Elevation 


F-< 


PRN  I 


SV  RISE/ST  GRAPH 


31” 

29T 

28T 

27T 

2T 

2S~~ 

24T 

23T 

22" 

21"' 

20" 

19"' 

18" 

ITT 

1ST 

1ST 

14" 

12"' 


4:30  5:00  5:30  6:00 


A 

4:00 


Figure  F  -  4.  SV  Rise/Set  and  Number  of  Visible  SVs 
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Figure  F  -  5.  Elevation/Time 
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Appendix  G 


Plots  of  North,  West  and  Azimuth  Tilts  (Figure  G  -  1  shows  conservative  Q-tuning, 
Figure  G  -  2  shows  non-conservative  Q-tuning);  Plots  of  North,  West  and  Vertical 
Velocity  (Figure  G  -  3  shows  conservative  Q-tuning,  Figure  G  -  4  shows  non¬ 
conservative  Q-tuning). 


Plot  Legend: 

••  true  error  (mean  error  ±  ctrue) 
filter  predicted  error  (0  ±  Ofjjter) 
mean  error 
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Figure  G  -  2.  North,  West  and 


rime  (sec) 

Jth  Tilts  (Non-Conservative  Q-tuning) 


1-3 


Vert.  Velocity  Error  (ft/sec)  West  Velocity  Error  (ft/sec)  North  Velocity  Error  (ft/sec) 


1 1 - 1 - i - 1 - 1 - 1 - 1 - r 


-|  I _ i _ i _ i _ i _ i _ i _ i _ 

0  500  1000  1500  2000  2500  3000  3500 


G-4 


West  Velocity  Error  (ft/sec)  North  Velocity  Error  (ft/sec) 


- 1 - 

— 

- ! - 

- 1 - 

/ — 

**  .  .  •  • 

r: . . 

■  _ 

- -  — — 

•-> 

. 7 

_  _  _  _  ^  •  ■*  _  * 

_ 1 - 

_ 

_ 1 _ 

- 1 _ 

■j  i _ i - 1 - 1 _ i _ i _ i _ i _ 

0  500  1000  1500  2000  2500  3000  3500 


. _ .  2 1  i  i  i  i  i - 1 - r 

o 


,2i _ i - 1 - 1 _ i _ i _ i _ i_ 

0  500  1000  1500  2000  2500  3000  3500 

Time  (sec) 


Figure  G  -  4  North,  West  and  Vertical  Velocity  (Non-Conservative  Q-tuning) 
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Appendix  H 


Plots  of  Case  I:  Barometric  Altimeter,  0.4  nm/hr  INS,  P-Code  GPS  Using  the  Tanker 
Flight  Profile. 

Plot  Legend: 

. true  error  (mean  error  ±  atrue) 

filter  predicted  error  (0  ±  afjjter) 
mean  error 


Error 
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Figure  H  -  1.  Latitude  and  Longitude  Error 
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Figure  H  -  2.  Aircraft  Altitude  and  Baro-Altimeter  Error 
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Figure  H  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  H  -  5  North,  West  and  Vertical  Velocity  Errors 
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Appendix  I 


Plots  of  Case  II:  Barometric  Altimeter,  0.4  nm/hr  INS,  P-Code  GPS  and  Radar  Altimeter 
Using  the  Tanker  Flight  Profile. 

Plot  Legend: 

. true  error  (mean  error  ±  atrue) 

filter  predicted  error  (0  ±  ctfiiter) 
mean  error 
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Figure  I  -  1.  Latitude  and  Longitude  Error 
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Figure  1-4  North,  West  and  Azimuth  Tilt  Errors 
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Figure  1-5  North,  West  and  Vertical  Velocity  Errors 


1-6 


GPS  User  Clock  Drift  (ft)  GPS  User  Clock  Bias  (ft) 


Time  (sec) 


Figure  1-6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Appendix  J 


Plots  of  Case  III:  Barometric  Altimeter,  2.0  nm/hr  INS,  P-Code  GPS  Using  the  Tanker 
Right  Profile. 

Plot  Legend: 
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Figure  J  -  1.  Latitude  and  Longitude  Error 
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Figure  J  -  2.  Aircraft  Altitude  and  Baro-Altimeter  Error 
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Figure  J  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  J  -  5  North,  West  and  Vertical  Velocity  Errors 
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Figure  J  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Appendix  K 


Plots  of  Case  IV:  Barometric  Altimeter,  2.0  nm/hr  INS,  P-Code  GPS  and  Radar 
Altimeter  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  K  -  2.  Aircraft  Altitude  and  Baro- Altimeter  Error 
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Figure  K  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  K  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Plots  of  Case  V:  Barometric  Altimeter,  4.0  nm/hr  INS,  P-Code  GPS  Using  the  Tanker 
Flight  Profile. 

Plot  Legend: 
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Figure  L  - 1.  Latitude  and  Longitude  Error 
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Figure  L  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Appendix  M 


Plots  of  Case  VI:  Barometric  Altimeter,  4.0  nm/hr  INS,  P-Code  GPS  and  Radar 
Altimeter  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  M  - 1.  Latitude  and  Longitude  Error 
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Figure  M  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 


M-4 


Error  (arcseconds) 


M-5 


Vert.  Velocity  Error  (ft/sec)  West  Velocity  Error  (ft/sec)  North  Velocity  Error  (ft/sec) 


M-6 


GPS  User  Clock  Drift  (ft) 


Time  (sec) 

Figure  M  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Plots  of  Case  VII:  Barometric  Altimeter,  0.4  nm/hr  INS,  P-Code  GPS  and  Single 
Pseudolite  Using  the  Tanker  Flight  Profile. 
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Figure  N  - 1.  Latitude  and  Longitude  Error 
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Figure  N  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  N  -  4  North,  West  and  Azimuth  Tilt  Errors 
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Figure  N  -  5  North,  West  and  Vertical  Velocity  Errors 
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Figure  N  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Appendix  O 


Plots  of  Case  VD3:  Barometric  Altimeter,  2.0  nm/hr  INS,  P-Code  GPS  and  Single 
Pseudolite  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  0-1.  Latitude  and  Longitude  Error 
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Figure  0-3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Appendix  P 


Plots  of  Case  XI:  Barometric  Altimeter,  4.0  nm/hr  INS,  P-Code  GPS  and  Single 
Pseudolite  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  P  -  2.  Aircraft  Altitude  and  Baro-Altimeter  Error 
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Figure  P  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  P  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Plots  of  Case  X.  Barometric  Altimeter,  0.4  nm/hr  INS,  P-Code  GPS,  Single  Pseudolite 
and  Radar  Altimeter  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  Q  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  Q  -  4  North,  West  and  Azimuth  Tilt  Errors 
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Plots  of  Case  XI:  Barometric  Altimeter,  2.0  nm/hr  INS,  P-Code  GPS,  Single  Pseudolite 
and  Radar  Altimeter  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 
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Figure  R  - 1 
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Figure  R  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Figure  R  -  5  North,  West  and  Vertical  Velocity  Errors 
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Figure  R  -  6  GPS  User  Clock  Bias  and  GPS  User  Clock  Drift 
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Appendix  S 


Plots  of  Case  VII:  Barometric  Altimeter,  4.0  nm/hr  INS,  P-Code  GPS,  Single  Pseudolite 
and  Radar  Altimeter  Using  the  Tanker  Flight  Profile. 

Plot  Legend: 

. true  error  (mean  error  ±  cr^g) 

filter  predicted  error  (0  ±  Ofjjter) 
mean  error 


S-l 


0  500  1000  1500  2000  2500  3000  3500 


Time  (sec) 

Figure  S  -  1.  Latitude  and  Longitude  Error 


S-2 


Aircraft  Altitude  Error  (ft) 


Aircraft  Altitude  Error  (ft)  Longitude  Error  (ft)  Latitude  Error  (ft) 


20 1 - 1 - 1 - 1 - 1 _ i _ i _ 

3600  3650  3700  3750  3800  3850  3900 


Time  (sec) 

Figure  S  -  3  Latitude,  Longitude  and  Aircraft  Altitude  Error 
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Appendix  T 


Plots  of  Cases  XIII,  XIV  and  XV:  Barometric  Altimeter,  0.4  nm/hr,  2.0  nm/hr  and  4.0 
nm/hr  INS,  P-Code  GPS  Using  the  Tanker  Flight  Profile  with  A  Single  GPS  Outage. 
Note  that  all  three  filters  seem  to  believe  the  GPS  clock  bias  error  will  blow-up  during  an 
outage,  while  the  true  GPS  clock  bias  error  seems  to  stay  unchanged  due  to  it  being 
modeled  as  a  random  constant  bias  in  the  GPS  truth  model. 

INS  plot  order  for  Figures  T  -  1  to  T  -  12  is  as  follows: 

•  0.4  nm/hr  (top  plot) 

•  2.0  nm/hr  (middle  plot) 

•  4.0  nm/hr  (bottom  plot) 


Plot  Legend: 
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Figure  T  -  9  West  Velocity  Errors 
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Figure  T  -  10  Vertical  Velocity  Errors 


T-ll 


GPS  User  Clock  Bias  (ft)  GPS  User  Clock  Bias  (ft)  GPS  User  Clock  Bias  (ft) 


)'■ - 1 - ‘ - I - 1 - 1 - - - - I _ L_ 

0  500  1000  1500  2000  2500  3000  3500 


T-12 


GPS  User  Clock  Drift  (ft)  GPS  User  Clock  Drift  (ft)  GPS  User  Clock  Drift  (ft) 


M3 


Appendix  U 


Plots  of  Case  XVI:  Barometric  Altimeter,  4.0  nm/hr  INS,  P-Code  GPS,  Single  Pseudolite 
and  Radar  Altimeter  Using  the  Tanker  Flight  Profile  with  Two  GPS  Outages. 

Plot  Legend: 
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